APPLICATION OF A FINITE ELEMENT METHOD TO 
THE BAROTROPIC PRIMITIVE EQUATIONS 



Donald Ernest Hlnsitian 



> 2318 ***' 






NAVAL POSTCRADOATE SCHOOL 

Monterey, California 










TH 




APPLICATION OF A 
FINITE ELEIffiNT METHOD TO THE 
BAROTROPIC PRIMITIVE EQUATIONS 

by 

Donald Ernest Kinsman 



September 1975 

Thesis Advisors: G.J. Haltiner 

R.T. Williams 



Approved for public release; distribution unlimited. 



Tl704f)3 



UNCLASSIFIED 



SECUBITV CLASSIFICATION OF THIS PACE fl«>.n Oa>« Enltrtd) 



REPORT DOCUMENTATION PAGE 


READ INSTRUCTIONS 
BEFORE COMPLETING FORM 


1. REPORT NUMBER 


2. GOVT ACCESSION NO. 


3. RECIPIENT’S CATALOG NUMBER 


4. title Cand 5u6(///«J 

Application of a Finite Element Method 
to the Barotropic Primitive Equations 


5. TYPE OF REPORT 4 PERIOD COVERED 

Master's Thesis; 
September 1975 


8. PERFORMING ORC. REPORT NUMBER 


7. autmoR('«; 

Donald Ernest Kinsman 


8. contract or grant number^*; 


9. performing crcanization name And address 

Naval Postgraduate School 
Monterey, California 93940 


10. program element, projec*^. task 

AREA & WORK UNIT NUMBERS 


11. controlling office name and ADDRESS 

Naval Postgraduate School 
Monterey, California 93940 


12. REPORT DATE 

September 1975 


13. number of pages 
116 


U. monitoring agency name & ADDR ESS(’i/ di//ar*nl from ControUing 0//ica; 


IS. SECURITY CLASS, (ol Ihit ripon) 

Unclassified 


15«. DECLASSI FI cation/ down GRADING 
SCHEDULE 



16. distribution statement Co/ <hi» R»pofO 



Approved for public release; distribution unlimited. 



n. Distribution statement (ot th» »b»trmct Biock 20 , a dHi»r»nt from Report) 



18. SUPPLEMENTARY NOTES 



IS. KEY WORDS (Continu* on fvrf if nocooomry mnd Idontify by biock numbor) 

Barotropic Primitive Equations 
Finite Element Method 
Icosahedral Grid 



20. ABSTRACT (Continue on r*v«r*« •id* if nocmaaory mtd idontity by biock numbor) 

A finite element application to the barotropic primitive 
equations is presented including theoretical development and 
the model used. Analytic initial data is generated in order 
to verify as well as possible the accuracy of the model. A 
comparison of the model with similar finite difference schemes 
shows that this finite element method exhibits better phase 
speed propagation than comparable second and fourth order 



DD I JAN *73 1473 EDITION OF 1 NOV 66 1$ obsolete 
(Page 1) S/N 0103-0M-660I | 



UNCLASSIFIED 

SECURITY CLASSIFICATION OF THIS PACE 



UNCLASSIFIED 

^iiCvjv^Ty Classification of this 



(20. ABSTRACT Continued) 

finite differencing and is competitive in the size of the 
allowable time step. 



DD Form 1473 
, 1 Jan 73 

S/N 0102-014-G601 



UNCLASSIFIED 



2 



SEConiTV CLASSIFICATION OF THIS PAOEr»T<»'' 



Application of a 
Finite Element Method to the 
Barotropic Primitive Equations 

by 

Donald Ernest BJinsman 
Lieutenant, United States Navy 
B.S., United States Naval Academy, 1968 



Submitted in partial 
requirements for 



fulfillment of the 
the degree of 



MASTER OF SCIENCE IN I>1ETE0R0L0GY 



from the 



NAVAL POSTGRADUATE SCHOOL 
September 1975 






ABSTRACT 

A finite element application to the barotropic primitive 
equations is presented including theoretical development and 
the model used. Analytic initial data is generated in order 
to verify as well as possible the accuracy of the model. A 
comparison of the model with similar finite difference 
schemes shows that this finite element method exhibits better 
phase speed propagation than comparable second and fourth 
order finite differencing and is competitive in the size of 
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I 



INTRODUCTION 



With the advent of the computer, there have been 
tremendous advances in the field of nvunerical weather 
prediction. Although there are several methods available for 
solving the prediction equations , the main thrust of the 
development of the numerical models has been with the finite 
difference method. Parallel developments in the engineering 
fields have used other methods as well as finite difference 
schemes. One recent technique being used is the finite 
element method. The purpose of this paper was to develop a 
finite element application to a barotropic primitive equation 
model. The objective was to learn the characteristics of a 
finite element model and compare it with a similar finite 
difference model. Analytic initial conditions were used to 
insure the most accurate analysis possible and to simplify 
the comparisons. 
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II. FINITE ELEMENTS 



Although the finite element method has only recently 
emerged as an efficient means of solving differential 
equations, its beginnings can be traced loosely back to early 
times. The basic concept of the finite element method is 
that a solution can be accurately determined by a sum of 
simple, easily computable functions. This procedure is not 
new. Martin (1973) states that a Chinese engineer in 
A. D. 480 was able to determine upper and lower bounds on 
the value of u to seven digits. This was accomplished by 
accurately determining the area of a circle with slender 
inscribed and circumscribed polygons. From calculus of 
variations, the classical Rayleigh-Ritz method shows how to 
approximate a solution to certain problems with a linear 
combination of any set of linearly independent functions . 

The Rayleigh-Ritz method determines the weights that are 
associated with each function. For a great many years both 
engineers and mathematicians were hung up on the idea that 
each of these functions m.ust be essentially non-zero over the 
entire domain of the problem. For a large complex solution, 
the determination of the weights could be prohibitive. Then 
in the early 1950 's engineers applied the Rayleigh-Ritz method 
subdividing the entire domain into many smaller pieces or 
elements. Hence the anachronism, finite elements. In this 
respect, the spectral method and the finite element method 
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ismm. 





are quite similar. The spectral method is a combination 
of sines and cosines defined over the entire domain while 
the finite element method is a combination of low order 
polynomials/ each polynomial or basis function, being non- 
zero only over a small "finite element." The question is 
then how to choose the coefficients for the polynomial to 
best approximate the ansv;er to the equation. The classical 
Rayleigh-Ritz procedure for linear self-adjoint problems 
formulates the problem as the solution to a minimization 
of a positive definite functional. The coefficients then are 
chosen to minimize the error in using a finite number of terms. 
Galerkin proved that the same coefficients result if the 
problem is formulated as a self-adjoint linear differential 
equation and then the coefficients are chosen in the following 
manner. Insert the linear combination of N basis functions 
into the equation, multiply this resulting equation by N 
linearly independent test functions (usually the basis 
functions themselves) and integrate to get N equations for 
the N coefficients. (Observe that if the test and basis 
functions are the same, multiplying the basis and test functions 
involves the integral of the squares of the functions , and 
thus this procedure is related to the least square error — a 
fact which Galerkin proved.) 

One may suspect that this same procedure would also give 
a reasonable error for nonlinear, non-self adjoint problems. 

In summary, the Rayleigh-Ritz-Galerkin procedure for 
nonlinear, non-self adjoint equations involves subdividing 



15 



the domain into a set of elements, approximating the dependent 
variables by a linear combination of low order polynomials 
(having value zero except over the particular element) , 
inserting these approximations into the equation, multiplying 
the equation by a test function (also a low order polynomial 
having a one-to-one correspondence with the basis functions) 
whose purpose is to minimize the residual between the 
approximate solution and the actual solution, integrating 
over the entire domain, and finally solving the system of 
equations assembled during the integration for the solution. 
Since the basis functions are zero over almost all of the 
area over which integration takes place, this procedure 
produces a system of equations for which the matrix is very 
sparse. The Rayleigh-Ritz-Galerkin method has become the 
most popular finite element method and is used in this paper. 
Martin (1973), Zienkiewicz (1971), Strang (1973), Schultz 
(1973) , Norrie (1973) , and Aziz (1972) give in-depth 
theoretical descriptions of the finite element method. 

The first step in the Rayleigh-Ritz-Galerkin method 
requires the division of the domain into a set of elements. 

The element shape in this paper was chosen to be triangular. 
Section III contains a description of the subdivisions 
utilized. 

The next step requires representing the dependent variables 
as a linear combination of low order polynomials. In this 
paper the polynomials used for both the test and basis 
functions were linear in A and 6. A linear basis function 
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over a particular element can be visualized as a plane with 
value one at one of the vertices of the triangle, value zero 
at the other two vertices of the triangle and identically 
zero over the rest of the domain. Since any one element has 
three points , there are three planes used in approximating 
a function over each element. This can be seen in Figure 1. 
Any function, say u, can be represented over the particular 
element as a linear combination of the three planes shown 
in Figure 1 by the equation 

’^element ' “j '^j 

where represents the scalar value of u at point j of the 
element and are the basis functions . At the boundary of 
two adjacent elements, the dependent variables are continuous. 
It may be noted, however, that since only linear functions 
are used for the representation, derivatives are not 
continuous along the boundaries. 

After the approximation for the variables have been 
substituted into the equations, the next step involves 
multiplying by a test function. Since in this paper, the 
test and basis functions were the same. Figure 1 also shows 
the three planes making up the test functions over any one 
element . 

The next step in the Rayleigh-Ritz-Galerkin procedures 
is to integrate over the domain. Since each basis and test 
function is zero over the domain except over the particular 
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FIGURE 1. Three planes representing a fxinction over an 
element . 
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element, the global integration can be performed by 
integrating locally over each element. In this manner the 
equations are written and then solved for at each node. The 
basis and test functions may be visualized from the view- 
point of a node instead of an element as shown in Figure 2 
where the six basis functions having value one at point i,j 
for each of the six elements are shown. The six planes 
(basis functions) shown make up a "pyramid function" for 
the grid point i,j. This "pyramid" function is the test 
function that multiplies the variable representation as 
described above. If this procedure is repeated for the N 
grid points a system of equations with N equations and N 
unknowns will be generated. 



19 



HEIGHT = 1 







FIGURE 2. Pyramid function at point i,j 
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III. BAROTROPIC PRIMITIVE EQUATION MODEL 



The time integration was performed using a Galerkin finite 
element application to the so called shallow water, "primitive" 
equations. A time extrapolated Crank-Nicholson method was 
used during formulation. The Crank-Nicholson method is 
explained in some detail in Appendix A. 

A. PRIMITIVE EQUATIONS 

The primitive equations, in spherical coordinates, for 
this model are as follows ; 



8t ■ a ^ST-e W + a - co — - e W <♦'' 






8u ^ u 3u ^ V 8u 
3t a cos e 3A a 36 



— tan 6 - 2n sin 6 v 



a cos 0 3 A 



3 4> _ 



O, LO— ^ 

(III-2) 



3v , u 3v , V 3v , u . n _L • /N 

•s-r + K -rr + “ + — tan 6 + sin 6 u 

3t a cos 6 3A a 30 a 



+ i = 0 

a 36 



(III-3) 



After equations (III-l) , (III-2) , and (III-3) were 
expressed in finite element form, they were solved in the 
following order: the continuity equation (III-l) , the 

u-equation (III-2) , and the v-equation (III-3) . Each equation 
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was solved by an iteration procedure until a converged 
solution was achieved before proceeding to the next equation. 

B . GRIDS 

Two different grids were used for experiments during 
development of the computer model, namely, the simple X,B grid 
and the icosahedral grid. The A, 6 grid conserves angular 
distance in radians while the icosahedral grid nearly 
conserves linear distance in meters. 

1 . A, 9 Grid ' 

Since the primitive equations are normally expressed 
in the spherical coordinate form, it seemed natural to use 
a grid simply subdivided by longitude and latitude. 

. 7T 7T 

Longitude ranged from zero to 2 tt and latitude, from ~ 2 2 ’ 

The grid was then evenly si±»divided into the desired intervals . 
A ten-degree interval generated a "square" with a side of 
length ten degrees. The construction of one diagonal across 
each square further subdivided each square into two triangles 
with the same square radian area. (See Figure 3.) The 
intersection of adjacent triangles' apices defined a node. 

Each node was then assigned a global correspondence number. 

The nximbering scheme ran along constant latitude lines in 
order to keep the bandwidth of the coefficient matrix as 
narrow as possible. Cyclic continuity was maintained by 
numbering the nodes at A = 2tt, with the same number as the 
A = 0 node given the same latitude circle. Every internal 
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node was supported by six other nodes; the north and south 
pole nodes were supported by four other nodes. A drawback 
to this grid was the north and south poles being represented 
by a line rather than a node. The primitive equations in 
the spherical coordinate form are singular at the poles. 

Thus this grid placed many nodes at one point (the north or 
south pole) where the equations are singular. 

2 . Icosahedral Grid 

Williamson (1968) and Sadourny, et al. (1967) 
pointed out the advantages of a grid which is nearly 
homogeneous with respect to area, and both described methods 
of generating a grid which used an icosahedral spherical 
figure. Cullen (1974) also used a regular icosahedron while 
integrating the primitive equations . 

A regular icosahedron on a sphere consists of twenty 
major spherical triangles with twelve vertices (see Figure 4) 
An icosahedral grid is then superimposed by subdividing the 
major triangles into smaller triangles of nearly uniform 
area, which is the most important feature of this grid. If 
a model can be run on a global grid with nearly equal area 
triangles then the problem of decreasing time steps 
associated with decreasing Ax as the poles are approached 
can be eliminated. Cullen (1974) subdivided each major 
spherical triangle along latitude/longitude circles but 
pointed out that a great circle subdivision would best 
conserve equal areas. 
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FIGURE 4. Regular icosahedron 
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The icosahedral grid used by this model consisted 
of a subdivision of each of the twenty major spherical 
triangles. Each major spherical triangle side was divided 
into n segments and the points were joined by arcs of great 
circles to produce a grid as shown in Figure 5. The 
mathematics are shown in greater detail in Appendix B. The 
number of points, N, in the grid where each side of a major 
spherical triangle is divided into n pieces is given by the 
following formula 

N = lOn^ + 2 (III-4) 

The number of triangles, T, generated by the same grid is 
given by 

T = 20n^ (III-5) 

With the exception of the vertices of the major spherical 
triangles, each node was supported by six other nodes. At the 
twelve major vertices, the nodes were supported by five nodes. 
The global correspondence number given each node was assigned 
by starting. at the north pole and going around the latitude 
band as shovm in Figure 6. The points of the triangles lying 
on zero and 360 degrees longitude were assigned the same 
global correspondence number to maintain cyclic continuity. 

This can be seen in Figure 7. 

3 . Global Correspondence Table 

A Galerkin finite element approach requires inner 
products of the dependent variables, say f, with a pyramid 
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FIGURE 5. Icosahedral grid 
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FIGURE 6. Global correspondence numbering scheme 
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FIGURE 7. Elements displaying cyclic continuity 
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function, . The inner product in spherical coordinates 
is defined as 

<f(X,0),V. > = // fU,6)V. a^ cos 0 dAd0 (III-6) 

global 

The pyramid function at any point i, has value one at i, 

linearly decreases to zero at the surrounding points and 
remains zero over the remainder of the globe. The global 
integration can be performed by integrating the particular 
function on each triangle, and placing the value of the 
integral in the proper place in the global matrix. A further 
explanation of the interaction between basis functions and 
test functions can be found in Section III.C.l. The global 
correspondence table is a means of properly scattering each 
triangle's surface integral into the appropriate space in 
the global coefficient matrix. 

Given a triangle with local coordinates 1, 2, and 3, 
the inner product over the triangle will produce a value for 
the three points 1, 2, and 3. Point 1 will receive the 
values from the interplay with itself (1,1), with point 2 
(1,2) and with point 3 (1,3) . Point 2 will receive values 
from interplay with (2,1), (2,2) and (2,3). Similarly point 
3 will receive values from (3,1), (3,2) and (3,3). Each 

triangle has a local 3x3 matrix which needs to be scattered 
into its global position in the coefficient matrix. Since 
each of the N points are multiplied by N global pyramid 
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functions, an N x N global coefficient matrix is formed 
during a complete integration over the globe for all points. 
Each row, i, in the matrix represents the equation for point 
i. A table had to be set up which correlates the local 3x3 
matrix to its place in the global N x N matrix. The first 
step in developing a global correspondence table is to number 
all the nodes. The numbering scheme should be chosen so as 
to minimize the bandwidth of the coefficient matrix. For 
this paper, another approach was taken which reduced the 
storage requirement to a minimum (Section III.C.3). George 
(1971) includes a siibroutine which, using any grid, will 
optimumly arrange the correspondence table to provide minimvira 
bandwidth . 

The icosahedral grid is indexed by starting at the 
north pole, point 1, and numbering around the next consecutive 
latitude band until all latitude bands are completed. The 
next step is to impose a local 1-2-3 coordinate onto each 
triangle and compile a local versus global table. Table I 
shows the first 10 triangles and its global correspondence 
table. The correspondence table is kept in a matrix and 
called during area integrations to scatter the 3x3 local 
matrix into the proper place in the coefficient matrix. The 
Computer Program section contains a program to generate the 
correspondence table. 
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10 


10 


3 


10 


4 



TABLE I. Global correspondence nui'nbers for 
the first ten triangles 
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C. FINITE ELEMENT APPLICATION 



K 

£, 



1 . Area Integrations 

The inner product <Vj,V^> can be computed by performing 
the necessary integration over each triangle separately and 
distributing the values to the proper position in the matrix. 

The actual integration can be done by a numerical scheme such as 
Gauss quadrature, or by use of an analytic approach. Due to 
the low order polynomials (linear) used, exact expressions 
may be derived for the integration of any function over an 
element. These formulas are used through a coordinate system 
called area coordinates. Strang and Fix (1973) state that 
area coordinates are known to engineers as triangular 
coordinates and to mathematicians as barycentric coordinates. 
Desai and Abel (1972) call them natural coordinates. 

In the formulation, it becomes necessary to perform 
the inner product <Vj,V^>. Given the three pairs of 
coordinates (x,y), of the triangle in the x,y plane, there are 
three corresponding coordinates (Lj^,L 2 ,L 2 ) in the natural 
coordinate system. This is shown in Figure 8. Zienkiewicz 
(1971) shows the relationship between a point x,y and the 
area coordinates (Lj^^,L 2 ,L 2 ) , as 

X = Lj^Xj^ + 1^2X2 + ^3^3 (III-7) 

y = + 1*2^2 ^3^3 (III-8) 

1 = Lj^ + L 2 + L^ (III-9) 

where L^^ = h^/A , L 2 = A 2 /A , L^ = A^/A 
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and A = total area of triangle; and A^,A 2 »A 2 = area of the 
smaller triangles. He also shows that 

// dxdy = — 2A . (III-IO) 

A ^ (a+b + c + 2)l 



This formula was used to perform the integration. For 
example, the inner product <Vj,V^> at j=i=l is 



<V, ,V^> = // V ^ dxdy 
.A 



21 0 ! 01 
(2 + 0 + 0 + 2 ) 1 



2 

24 



2A 



A 

6 



and the inner product 



<Vj,V.> 



at point j=2, i=l is 



<V 



2 



^1> 



// V V dxdy 
A ^ 



1 ! 11 01 
(1 + 1 + 0 + 2 ) 1 



2A 



1 OA A 

‘ 24 ' “ 12 



With this formula, any necessary product of basis 
functions can be integrated over one triangle. 

If equations (III-7) , (III-8) , and (III-9) are now 

solved foi L^, L 2 , and written in matrix form the result 
is 




(III-ll) 
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FIGURE 8. Cartesian coordinates vs. natural coordinates 
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FIGURE 9 . Triangle definitions for area coordinates 
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as shown in Desai and Abel (1972) where the a's and b's are as 
defined in Figure 9. Differentiation of (III-ll) shows that 



3X 




b. 

1 

2A 




(III-12) 



1_ 

ay 






(III-13) 



These two equations will be used when evaluating the 
derivatives of the basis and test functions. 



For example, assuming V. = L. , the derivative V. 

9 9 9x 



at point 1 is 



V 



bi 3V^ av^ b^ 

lx "" 2A aL^ 7K 2A 9 L^ 



but 



3V, 



= 0 



aL. 



= 0 



av 

and ’5^“ ~ ^ point 1, 



Consequently , 



^Ix 2A 



The inner product point j=2, i=l becomes 






A 

b. 



11 0! 0! 



7a (1 + 0 + 0 + 2 ) ! 
^2 

T* * 
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Thus with these simple formulas, all the inner products can 
readily be evaluated. 

2 . Nonlinear Terms 

The nonlinear terms in the system of equations were 
quasi-linearized by the time extrapolated Crank-Nicholson 
method. VThen a uu^^ term was encountered, it was replaced 
by u*u, where 



u* 




- i 



(III-14) 



Thus when solving for time level N+1, all the * quantities 

will be known. In this manner the nonlinear terms are quasi- 

lienarized in time to facilitate integration. The inner 
* 

product <a . V . . , V. > represents u*u, in the notation of 
this paper. In matrix notation this becomes 



a.A. . = a . <a*V, V. , ,V. > (III-15) 

3 13 3 k k 3X' i 

In the local 3x3 matrix, any point within the local matrix 
becomes 



Ai- = 









(III-16) 



In this manner the known * functions are integrated into the 
Galerkin space. The nonlinear terms are "linearized." Each 
equation becomes one equation v/ith one unknown. The three 
equations are coupled when the new information at time level 
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N+1 is used to solve the equations at time level N+2. A 
constraint on a linearization scheme is that a large time 
step cannot be taken. The time step should not be so large 
as to cause the change in the variable to be larger than the 
truncation error of the approximation. This linearization, 
which uncouples the three equations, causes some inconsistency 
with respect to the three dependent variables. 

3 . Matrix Storage 

During a global integration, an N x N coefficient 
matrix is generated which is very sparse due to the fact that 
the maximum number of triangles supporting any one point is 
six. These six triangles provide interaction between the 
seven points involved. Each row, i, in the N x N matrix 
therefore represents the equations written down at point i, 
and each row would have, at a maximum, seven entries. If 
the matrix were condensed by omitting all terms identically 
zero, the size would be N x 7. This represents a sizable 
reduction in core storage. 

The method of condensed matrix storage offers the 
advantage of minimum core storage. For large fields, a 
complete N x N matrix can easily exceed core capabilities 
of most computers. The storage of a sparse matrix would be 
extremely wasteful, hence, the condensed method was adopted. 

If the coefficient matrix were symmetric or had a narrow 
bandv/idth, then storage could be easily accomplished using 
a smaller amount of computer core. The icosahedral grid 
offers neither a symmetric matrix nor one with a narrow 
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bandwidth. The disadvantage of the matrix storage scheme 
was the need to search a correlation matrix when a value 
needed to be placed in the coefficient matrix. An efficient 
searching routine was developed and can be seen in the 
Computer Program section under Subroutine SEARCH. But, 
unfortunately, this subroutine had to be executed many times 
during the time integrations. 

In this model, an iterative procedure was used to 
solve the coefficient matrix. In order to reduce the N x N 
matrix to an N x 7 matrix, a correlation table had to be 

developed. The correlation table, a separate N x 7 matrix, 

was assembled which contained the seven points involved in 
any one row of the coefficient matrix. An example will show 
more clearly the process involved. Table II and Figure 10 
sho\^7 the triangle, points and correlation table involved with 
point 2. In an N X N matrix, point 2's equation (row 2) would 
have non-zero entries in position (2,1), (2,2), (2,3), (2,6), 

(2,7), (2,8), and (2,16). For example, in the N x 7 coeffi- 
cient matrix, entry (2,16) would be stored in (2,7) of the 

condensed matrix and entry (2,8) would be stored in (2,6) . 
During matrix multiplication, the proper address of a vector 
component had to be selected as well as the proper address 
in the N x 7 matrix. This can be accomplished by calling, 
for example, position (2,7) of the correlation matrix to 
find that entry (2,7) of the coefficient matrix goes with 
entry 16 of the vector. The process is reversed when it 
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Legend 



0 Triangle number 1 
1 Point number 1 




O’ 

Figure 10. Diagram of triangles 1, 5, 6, 7, 19, 20 and 
nodes 1, 2, 3, 6, 7, 8, 16. 



TABLE II. Correlation table for nodes 2, 7, 8 
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becomes necessary to scatter the local 3x3 matrix built 
during each triangle's surface integration. The global 
number of each point in the triangle determines which row 
into the N x 7 coefficient matrix the value will go. The 
colxomn is found by doing a search of the row from the 
correlation matrix until the second number is found. The 
colvunn of the correlation table in which the second number 
was found determines the column in the N x 7 coefficient 
matrix. For example, triangle 6 will have a local 3x3 
matrix with global numbers which looks like 



(2,2) 


(2,7) 


(2,8) 


(7,2) 


(7,7) 


(7,8) 


(8,2) 


(8,7) 


(8,8) 



Row 2 of the coefficient matrix will contain the values of 
(2,2), (2,7), and (2,8). A search of correlation matrix's 

row 2 shows that (2,8) of the N x N matrix goes into (2,6) of 
the N X 7 matrix. 

The correlation matrix is developed by searching the 
global correspondence table for the six triangles that contain 
point i. These six triangles are then compared and sorted to 
produce the seven points interacting with point i. The 
seven points are then arranged in ascending order. 

A copy of the program which produced the correlation 
table can be found in the Computer Program section. 
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IV. INITIAL CONDITIONS 



The initial conditions used in this paper were computed 
from an analytic solution to the nonlinear balance equation. 

Haurwitz (1940) developed a stream function and showed 
that harmonic waves computed from this stream function will 
move with constant angular velocity in a non-divergent baro- 
tropic atmosphere. The stream function, ip, used by Haurwitz 
is given by 

2 

ip = A* sin(mX - vt) sin 9 cos m9 - Ba sin 6 (IV-1) 

where A* and B are constants, m is the wave number, v is the 
angular wave velocity, a is the radius of the earth. 

The constant B is related to the wave number and angular 
wave velocity by 

V ^ B M(M+1) - 2 _ 2Q . 

m M(M+1) M(M+1) ' ^ 

where (v/m) is the angular phase speed, Q is the angular 
velocity of the earth and M = m+1 . To allow a phase speed 
propagation comparison with a finite difference scheme, the 
constant A* was arbitrarily chosen to be the same as the A 
in a M.S. thesis by Monaco (1975). Phillips (1959) used the 
Haurwitz stream function, ip, in the nonlinear balance equation 
to determine (}>. The nonlinear balance equation is 
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(IV-3) 



= fS7^4> - Vip-Vf - 2(|^)2 + 2 (^ 

8x^ 3y2 



which in spherical coordinates becomes 



1 r 1 9 <}) , 1 9 , ^ 9(j), , 

-y [ 5 7 + (cos e -^) ] 

a cos 6 8X cos 0 30 ° 



“T t 1 — ^ — (cos 0 ||) ] 

a cos^ 0 9X cos 0 30 



+ J- li M _ 2( i 

a^ 30 30 a^ cos 0 3A30 



ijL.)2 



2 _9?1 

a^ cos^ 0 3X^ 30^ 



(IV-4) 



Phillips (1959) found the solution to the nonlinear balance 
equation for the geopotential distribution, (p, to be 



<j) = a^ A(0) + a^ B(0) sin(mX) + a^C(0)(2 sin^ mX - 1) 



(IV-5) 



where 

A(0) 



2m 



|(2fi+B)cos^0 + i(^)^cos 0 [(m+l)cos^0 + (2m^-m-2) 



2m 



cos 0 



-] 



(IV-5.1) 



A* 

2(fi+B) ~ 

B(0) = (nt+i) (m+ ' 2 ' ) ~ cos"'0 ((m^+2m+2) - (m+1) ^cos^ 0] 



(IV-5. 2) 
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( 



and 



C{6) 



i(— 

4' 2’ 

a 



cos^^ e [(m+l)cos^0 - (m+2)] (IV-5.3) 



The foregoing geopotential disturbance, <p, must be added to 
the mean height for the level desired. The mean height was 
obtained from the NACA standard atmosphere (Haltiner and Martin, 
1957) . Although an attempt has been made to apply the model 
at the level of non-divergence, the equations used in the model 
permit divergence. Phillips (1959) found that the presence of 
divergence will cause the waves in the height fields to move 
slightly slower than a non-divergent model, especially for 
small wave numbers. 

Analytic initial winds can be obtained from the stream 
function, ip, as follows 



and 



u 



1 

a TF 



= 1 11 
a cos 6 dX 



(IV-6) 

(IV-7) 



From ip, in IV- 1, the values of u and v are 

u = - — [A*sin(mX-vt) cos^^6 - mA*sin(mX-vt) cos'^^S sin^G - Ba^cos 6] 

(IV-8) 

V = — [A*m sin 0 cos'^^0 cos(mX-vt)] (IV- 9) 

a 



These analytically determined winds from the nonlinear 
balance equation were chosen because they are in approximate 
gradient balance and hence more realistic than simply 
geos trophic winds. Subroutine SOLUT of the main program. 
Computer Program section, produced the initial fields for 
the latitude and longitude of each point. 
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V. WAVE ANALYSIS METHOD 



A harmonic analysis was performed around constant 
latitude circles to determine the phase speed and amplitude 
of meteorological waves. The initial fields were analyzed 
first to provide a reference condition. 

The Fourier series used was 

F(x) = A + y (A cos mx + B sin mx) (V-1) 

o m m 

m 

which can be represented as one trigonometric function given 
by 



F(x) = C + y C cos(mx-6 ) 
o i; m m 

m 



where 



m 



B 

m 

sin 

m 



A 



m 



cos 6_ 
m 



(V-2) 



(V-3) 



and 

= tan"^ ^ (V-4) 

In order to analyze phase speed propagation, nonlinear 
instability and energy conservation, numerous wave numbers, 
amplitudes and phase speeds were examined. 
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VI. GRID CONVERSIONS 



The grid that is most compatible to this finite element 
model is the icosahedral grid. But this grid is incompatible 
with the grids needed to perform harmonic analysis or 
graphical display. This fact could be a serious drawback to 
an operational finite element model. An advantageous 
characteristic of the finite element method is that the 
solutions are continuous over the element. For this model, 
linear fxinctions were used. To convert the icosahedral 
solutions to solutions at different points involves evaluating 
the function over the element at the point desired. This is 
where the finite element method will have an advcintage over 
the finite difference schemes. Finite difference schemes 
must interpolate from discrete point values to an interior 
point. In this finite element method, the basis functions 
are of the form 

f = d + bX + ee (VI-1) 

The value of the funciton is known at three points. These 
three equations can be written for points 1, 2, and 3 



*1 = “ 


+ 


bX^ 


+ 


e0j^ 


(VI-2) 




+ 


bX2 


+ 


e02 


(VI-3) 


^3 = ^ 
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bX3 


+ 


e02 


(VI-4) 
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We have a system of three equations and three unknowns. The 
value of the coefficients d, b, and e over ciny one element 
can be easily found. Once the coefficients are known, the 
value of the function at any point X,0 in the element can be 
found by evaluating the polynomial given in equation (VI- 1) . 

For this model a 72 x 35 regular 5® longitude- latitude 
grid was required for harmonic analysis and graphical output. 
A subroutine was developed which correlated each point in the 
72 X 35 grid to its appropriate triangle in the icosahedral 
grid. This is subroutine SURVEY in the main program in the 
Computer Program section. The value of the coefficients d, 
b, and e in equation (VI-1) were found for each triangle by 
subroutine EVAL in the main program in the Computer Program 
section. The value of the function at any point in the 
72 X 35 grid was evaluated by subroutine DISPL in the main 
program in the Computer Program section. 
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VII. EXPERIMENTS 



Experiment 1 . The model was run with the same conditions 
using the icosahedral grid and the X,9 grid. The purpose 
of this experiment was to compare a varying Ax length grid, 
as most non projection finite difference models use, with a 
near constant Ax grid such as the icosahedral grid. The 
conditions used were a 10-minute time step, wave number 4, 
amplitude times wave number = 28. x 10 meters and a phase 
speed of 10® longitude/day , 

Experiment 2 . The model was run on the icosahedral grid 
with wave numbers 4, 8, and 12. The phase speeds were all 
10® longitude/day. A 10-minute time step was used for all 
runs. The amplitude of the waves was changed so that the 
product of amplitude and wave number remained 28 x 10^ meters. 
This will constrain the order of magnitude of the north-south 
component to remain the same. The purpose of this experiment 
was to observe phase speed propagation of the various waves 
and compare them to the phase speed propagation characteris- 
tics of various finite difference schemes. 

Experiment 3 . The model was run with increasing time 
steps with wave number 4. The other conditions were the sam.e 
as in experiment number 1. Time steps of 10, 12, 15, and 
18 minutes were run for a full 72 hour prognosis. 
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VIII. RESULTS 



Experiment 1 . Initially the model ran on the A, 6 grid 
but after twelve hours the harmonic analysis near the polar 
regions showed increasing amplitudes in the high wave numbers. 
Eventually the solutions in the polar regions caused the 
model to become unstable. There are several possible 
explanations for this instability. From a stability criterion 
viewpoint, the Ax near the poles was so small using the A, 6 
grid that the 10-minute time step exceeded the stability cut 
off point. This is another major drawback of this grid. 

The decreasing Ax as the pole is approached is a common source 
of problems for normal finite difference schemes. From 
theoretical considerations, the finite element method should 
be well suited to an icosahedral grid and this should 
alleviate the converging Ax problem. The A,0 grid has 36 
points around each latitude circle. It can therefore resolve 
wave number 18. Any truncation error created while solving 
the matrix equations will appear as high frequency noise. 

Another reason is that the equations are singular at the 
poles. In order to avoid these singularities, the model does 
not predict u, v or <{> at the poles. Instead u and v are set 
to zero initially and are kept zero throughout the entire 
integration. A value for ({> at the poles is found by averaging 
the first latitude circle down from the pole and then setting 
the pole and the next latitude circle equal to the average 
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value. This will flatten the cap over the pole and make 
the u and v assumptions at the poles more consistent, although 
obviously this treatment of the poles is not completely 
desirable . 

Experiment 2 . Figures 11, 12, and 13 show the results 
of the phase speed propagation test for v/ave numbers 4, 8, 12 
respectively. While identical tests using finite difference 
schemes were not available, Maher (1974), in a Master's Thesis, 
performed phase speed propagation analysis with similar wave 
numbers and amplitudes using second and fourth order finite 
difference schemes with a compareible number of points . The 
comparison showed that the finite element method was more 
accurate than second or fourth order differencing for the 
cases examined. 

Due to the coefficient chosen, the amplitude was minimal 
in the waves at higher latitude. On the right side of each 
figure is a value showing the number of points around each 
latitude band. As the number of points per each latitude 
band decreases, the phase speed propagation decreases. The 
finite element method therefore has the same relationship as 
finite differencing where phase propagation is concerned. It 
takes more points per wave to improve the phase propagation. 
Charts A-I show the wave propagation in 12 hour increments 
for the three wave numbers tested. 

Experiment 3 . With the three equations being treated 
separately at each time level and the main forcing functions, 
Coriolis force, and pressure gradient force, being treated 
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explicitly, it is reasonable to expect a linear stability 
criterion approximately of the form 




where c = the phase speed of the fastest waves and Ax the 
minimum Ax in the grid. The actual stability criterion for 
this model on the icosahedral grid has not rigorously been 
worked out. However, in notes by Archer (1975) , the stability 
criterion for a simple wave equation was worked out and found 
to be comparable to that with similar finite difference 
schemes. The assumed stability criterion predicted a maximum 
time step of 19 minutes. It was found that the model with 
an 18 minute time step became unstable after 12 hours and 
after 24 hours with a 15 minute time step. However the model 
remained stable with a 12 minute time step up to 48 hours. 

In all cases, instability in the equatorial regions was 
observed after a sufficient period of integration. For the 
wave number 4 case, this instability occurred after 84 hours 
of real time integration. The instability occurred sooner 
for larger time steps. 
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LATITUDE 
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PHASE ANGLE (DEGREES LONGITUDE) 



FIGURE 11. Phase angle (degrees longitude) vs. latitude 
for icosahedral grid, wave nuinber 4, phase 
speed 10 “/day and A* = 7.0 x 10'. (Latitudes 
with near zero wave amplitude are not included 
and time is given in hours.) 
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NUMBER OF POINTS PER LATITUDE 



I 




Chart A. Initial PHI field analysis, wave nun±ier 4, phase 
speed 10 ’’/day, A* = 7.0 x 10^. 
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Chart B . 



12-hour PHI field forecast, wave number 4, phase 
speed 10°/day, A* = 7.0 x 10^. 
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Chart C. 



24-hour PHI field forecast, wave number 4, phase 
speed 10°/day, A* = 7.0 x 10^. 
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FIGURE 12. Phase angle (degrees longitude) vs. latitude for 
icosahedral grid, wave number 8, phase speed 
10°/day and A* = 3.5 x 10^. (Latitudes with 
near zero wave amplitude are not included and 
time is given in hours.) 
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NUMBER OF POINTS PER LATITUDE 




Chart D. Initial PHI field analysis, wave number 8, phase 
speed 10°/day, A* = 3.5 x 10'. 
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Chart E. 12-hour PHI field forecast, wave number 8, phase 
speed 10°/day, A* = 3.5 x 10^. 



59 




Chart F. 



24-hour PHI field forecast, wave number 8, phase 
speed 10®/day, A* = 3.5 x 10^. 
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13. Phase angle (degrees longitude) vs. latitude 
for icosahedral grid, wave number 12, phase 
speed 10°/day and A* = 2.3 x 10^. (Latitudes 
with near zero wave amplitude are not included 
and time is given in hours.) 
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Chart G. Initial PHI field analysis, wave number 12, phase 
speed 10°/day, A* = 2.3 x 10'. 
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Chart H . 



12-hour PHI field forecast, wave number 12, phase 
speed 10 “/day. A* = 2.3 x 10^. 
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Chart I. 24-hour PHI field forecast, wave nuiBber 12, phase 
speed 10°/day, A* = 2.3 x 10^. 
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IX. CONCLUSIONS 



This model showed that the finite element method is 
applicable to a meteorological system and is competitive with 
finite difference schemes. The size of the time step that 
can be taken by a finite element method is comparable with 
the one for a finite difference scheme. The accuracy of the 
phase propagation of the finite element method was closer to 
the analytic solution than second and fourth order finite 
differencing for the comparisons made. 

Based on the increased accuracy and comparable time steps, 
future research should be continued on the finite element 
method. Specifically, more advanced methods of solving the 
system of equations are available and should be utilized, for 
example, ADI or a cyclic reduction method. Also higher order 
polynomials should be used for the basis functions. In fact, 
the real advantage of finite elements is not realized until 
higher order polynomials are used, from which greater accuracy 
can be expected. 

Research should continue in the polar regions to ascertain 
the ability of the method to resolve flow over the poles. 

This will clearly improve the model as it will have realistic 
physical conditions at the poles instead of artificial 
boundary conditions. To extend the allowable time step, 
the equations should be written so that the three equations 
v;ould be coupled at any one time step. Also the instability 
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observed in the equatorial regions after long time integra- 
tions should be further investigated. 

Ihe barotropic model should be expanded into a multi- 
level baroclinic model with a heating package. Real data 
should also be applied to the barotropic and baroclinic 
models to test the method's ability to handle actual 
conditions . 

Finally, research should be performed using smaller 
elements perhaps selectively in certain areas. With the 
success with the icosahedral grid, the application of the 
finite element method to fine meshed models shows promise 
to provide greater resolution in specific regions. 
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APPENDIX A 



EQUATION FORMULATION 

In the initial stages of this paper several decisions 
were made as to the order and method of solving the three 
coupled equations. It was decided to uncouple each equation 
so that it would be one equation with one unknown. Each 
equation would be written in the time extrapolated Crank- 
Nicholson method-. Although this method has a slov;er 
theoretical convergence rate than solving the three equations 
coupled together, the number of computations per time step is 
reduced. If the three equations are coupled, then each 
equation's solution after one iteration at one time step is 
used in solving the other two equations for the same iteration 
for the same time step. This requires area integration of 
each equation for each iteration at any one time step. For 
example, with the Crank-Nicholson method the <j) in the u- 
equation is written at time N+J^ . Since this is an unknown, 
each iteration of the (|) equation for time, N+1, gives a 
different . The pressure gradient term in the u-equation 

would have to be integrated again as it is not the same as 
it was at the last iteration. The uncoupling of the three 
equations was chosen because it sharply reduces the number 
of computations per time step. The advantage of coupling the 
three equations at any one time step would be that the 
equations would be more accurate and consistent and larger 
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time steps would be possible. Even with the restriction on 
the size of the time step, a 12-minute time step was possible 
before instability was experienced. 

The first equation solved during one time step is the 
continuity equation 

t 

a cL e lx a ' cbs e le 



An inner product is now performed on each term with a 
global pyramid function . The inner product is defined 

<f(A/0)/V. > = jj f(A,0) V. a^ cos 0 dXd0 (A-2) 

global 

2 

The a is not carried since xt is a common factor and can be 
cancelled now. 

The continuity equation becomes 



94 > 

3t 



V. > 
1 



+ i 

a 



cos 6 3X + a "^cos 0 30 



^ (4>v cos 



0) /V^> 



= 0 



(A-3) 



If the second and third terms are integrated by parts, 
the continuity equation becomes 



< If ^ 



2tt , , ^ 3V. 

d0 — < — ^ , > 

a cos 0 3X 



+ — / ( 4>v cos 0 V.) 
a 1 



7T 
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3V. 
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IT a cos 0 3 0 



(A-4) 
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The second term in the above equation is zero due to 
cyclic continuity and the fourth term is zero due to the value 

TT TT 

of cos 6 at - j j . The continuity equation becomes 



<M 

at' 



V. > 



1 



i <— 

a cos 6 




1 ^ 4)V cos 6 ^^i ^ 

a cos 0 '36 



0 (A-5) 



To extrapolate the dependent variable forward in time the 
Crank-Nicholson method first uses a simple forward difference 
in time to N+1 followed by an average at time level N+1 with 
values at time level N for the space derivatives. The space 
derivatives, therefore, are evaluated at time level N+Jg. 

Since the three equations were uncoupled, any variable that 
is not the variable for the particular equation is evaluated 
at time N+% while the variable for the equation is evaluated 
at times N and N+1. The continuity equation becomes 
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(A- 6) 



The variables u and v at time N+% are represented by the 
second order approximation 



N+H * - 3 N 1 N-1 

u 2 = u*-ju - j u 



(A- 7) 



N+% _ * ~ 3 N 1 N-1 



(A-8) 
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This completes the uncoupling and the continuity equation 
becomes 



<(*N+1-4.N),V.> - ^ 

1 ZdL 
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(A-9) 



The continuity equation is now one equation in one unknown. 

Using Galerkin formulation with Einsteinian notation,^ 
the variables u, v and <{) are represented by 



u = 


0i.it) 


V. 

D 


(A-10) 


V = 


(t) 




(A-11) 


<i> = 


Yj(t) 


V. 

3 


(A- 12) 



where Vj are the low order polynomials that are the basis 
functions. In this paper we used linear polynomials of the 
form 



f = d + bX + e6 



(A-13) 



The coefficients a, 3 and y are the scalar values of the 
variables and are functions of time only. 



^A repeated index implies summation with respect to that 
index. 
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Using this new notation, the continuity equation becomes 






At 

2a 




a, V.V, 3V. 
k -| k 1 

cos 6 ' 3A 



N+1 * cos 6 

■^•=^3 H55-e 



At 

2a 



<Y 



N 



ajJVjV, 



. 3V. 

k 1. 



j cos 0 ' 3 A 



^ "^3®k 



N„* cos 6 



cos 0 



V.V, , 
D k' 



3V. 

X 

30 



= 0 



(A-14) 



where is the pyramid function at the point where the equation 
is being written. It is now clear why the pyramid function was 
introduced. This procedure is repeated for each grid point 
leading to a system of equations for the value of dependent 
variable at the N grid points. 

The trigonometric functions in the equation can be handled 
in several v;ays . Cullen (1974) treated the trigonometric 
functions as constants over each triangle. As the grid length 
decreases the accuracy of this approximation increases. The 
decision was made to interpolate the trigonometric functions 
into the Galerkin space. The trigonometric functions were 
represented by 

cos 0 = CjVj (A-15) 

sin 0 = HjVj (A-16) 
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For ease of notation the trigonometric functions were 
cancelled where possible from the inner products and v/e now 
define a new inner product 



2 tt/2 

<u,V. > = / / u V. dXde (A-17) 

^ 0 -it/2 ^ 



The continuity equation becomes 






N+1 * ^^i 



N+1 * ^^i 



2a 



N * ^^i 

<Y a, V.V, ,-rr^> 
— ' 3 k 3 k 3A — 



N * 



= 0 



(A-18) 



Since everything is known in this equation, with the exception 
N+1 

of Y • / and the Y- are functions of time only, we can take 

3 1 

the Yj outside the inner products and write the following matrix 
equation 



[A - At B] 



N 4 - n 1 



(A-19) 



where 

A 






and 



B = 



* 8V. * 3V. 

+ iekCLVjV3^V^,^> 



(A-19. 2) 
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The matrix equation is manipulated to solve for the 

u N 

change, e^ , in 









.V«,§|b-§|b) 



[A 

D 



2a 



B] = 



^ B1 + ^ B 

D J a 



e«(A - §|b1 



y" ^B 
D a 



where 



N+1 

j 






The zonal equation was developed in the same manner 
the continuity equation and the matrix equation derived 



Nr^ ^ Kl 

ej [A + D) 



At =T N At - 

— E - a . — D 

a Da 



where 



8V. 



^ = ^Vk8A 



3V. 






and 



* 






+ 



<a*p*nT V, V .v^ 
k 3 L k 3 L 



V. > 
1 — 



(A-20.1) 

(A-20.2) 

(A-20) 

(A-20.3) 

as 

was 

(A-21) 

(A-21.1) 

(A-21. 2) 
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The meridional equation becomes 





(A-22) 



where 




(A-22.1) 



and 





(A-22. 2) 



The matrices in equations (A-20) , (A-21 ) , and (A-22) 

are built by the procedures in Section III.C.l. Once the 
matrices for each equation have been built, the system of 
equations derived can be solved by any conventional process. 

A simple Gauss-Seidel iterative procedure was chosen to solve 
the matrix equations with a relative improvement check used 
as a cutoff criterion. Subroutine SOLVE in the Computer 
Program section solved the equations. 

While the purpose of this paper was to develop a finite 
element barotropic primitive equation model, it was hoped 
that a considerable savings in time could be realized. Further 
time savings can be achieved if more sophisticated techniques 
such as. Successive Over Relaxation, SOR, or Alternating 
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Direction Implicit, ADI, are used. Leslie and McAvaney (1973) 
show that the computing time can be greatly reduced by using 
some of the newer direct techniques. These latter were not 
employed during development since programs for these methods 
were not available in the library of programs. 
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APPENDIX B 



ICOSAHEDRAL GRID FORMULATION 



A spherical icosahedron is made up of twenty equilateral 
spherical triangles. Each interior angle is 2ir/5. The 
method of subdividing the icosahedron is accomplished in 
three steps. First, the top five major spherical triangles 
are partitioned. Second, the Southern Hemisphere triangles 
are reflected from the Northern Hemisphere. The last step is 
the subdivision of the equatorial triangles. 

The first Northern Hemisphere major spherical triangle 
is subdivided and the next four adjacent major spherical 
triangles are given the same solution for longitude and 
latitude of the nodes except that the longitudes are displaced 
the proper multiple of 2ir/5 radians. In all cases, the purpose 
of the subdivisions is twofold. The main purpose is to specify 
the latitude and longitude of every node. The second purpose 
is to assign a global correspondence number to each node. 

The second step in subdividing the sphere is to mirror 
the Northern Hemisphere spherical triangles into the Southern 
Hemisphere. The solutions are the same except that the 
longitudes are displaced tt/ 5 radians. 

The third step is the solution of the ten major spherical, 
equatorial triangles. The following description shows how 
one Northern Hemisphere triangle, two equatorial triangles 
and one Southern Hemisphere triangle can equally be subdivided 
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and solutions for longitude and latitude of each node be 
obtained (Figure 14). Three laws of spherical trigonometry 
are required. They are: 

Law of Sines 

sin a _ sin b 

sin A sin B ' 

Law of Cosines for Sides 

cos a = cos b cos c + sin b sin c cos A (B-2) 



Law of Cosines for Angles 

cos A = -cos B cos C + sin B sin C cos a (B-3) 

where a represents a side (spherical arc) of a triangle and 
A represents the corresponding opposite angle. 

Since all three angles of a triangle are known, the length 
of arc a can be determined by the law of cosines for sides. 
Since a is along a meridian, the length of this arc is the 
colatitude of point 2. The longitude of points 1 and 2 is 
arbitrarily chosen as zero radians. Arc a is then divided 
into n segments. The latitudes and longitudes of these points 
can easily be computed. Arc b can be segmented in the same 
manner. Arc d can be found by using the law of sines with 
angle C and arcs 3-1 and 4-1. Arc d is then segmented into 
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NORTH POLE 







t 



SOUTH POLE 



FIGURE 14. One Northern Hemispheric, two Equatorial and 
one Southern Hemispheric major spherical 
triangle. 
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the proper number of pieces. The three sides of triangle 
1-4-3 are knovm, therefore angle D can be found by the law 
of cosines for sides. Arc 1-3, arc e and angle D with the 
law of sines determine the colatitude of point 5, arc f. 

Arc 1-3, arc e and arc f with the law of cosines for sides 
determine angle E, point 5's longitude. All points in the 
Northern Hemisphere triangle can be solved for longitude 
and latitude by manipulating these laws, arcs and angles. 

All the points in a Southern Hemisphere triangle are 
mirror images of- those in a Northern Hemisphere triangle 
except that the longitudes are displaced it/5 radians and the 
latitudes have a minus sign. 

With angle (B+F) , arc 1-6, angle C/2, and the law of 
sines, arc g can be found. This arc is then divided into 
n segments. With angle (B+F), arc a, arc 2-7, and the law 
of cosines for sides, arc 1-7 can be found. This arc is the 
colatitude of point 7 but care must be taken as this arc 
can extend across the equator. Arc 1-7, angle (B+F) , 
arc (2-7) and the law of sines determine angle (2-1-7) , 
the longitude of point 7. Point 8's longitude and latitude 
can be foxind by a similar computation using the arcs and 
angles on the other side of the triangle. Arc 7-8 can now 
be found and subdivided. Angle (8-7-1) can also be found. 

Then triangle (9-7-1) can be solved to determine the latitude 
and longitude of point 9. The rest of the points in both 
major spherical equatorial triangles can be solved for latitude 
and longitude in a similar manner. The Computer Program Section 
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contains a program written in Fortran IV for an IBM 360 
computer which will subdivide an icosahedron simply by giving 
it the desired number of segments, n. 
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oooooooooooooooooooooonooooooooooooooooooooooooooooooooooooooo 



COMPUTER PROGRAMS 



ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

THIS IS THE MAIN BODY OF THE FINITE ELEMENT 
PROGRAM. 

THE INTEGER’i'Z STATEMENTS ARE USED TO STORE 
THE REQUIRED TABLES USED DURING INTEGRATION. 

SINCE THE LARGEST INTEGER STORED DID NOT 
EXCEED THE CAPABILITIES OF THE INTEGER*2 MODE 
IT WAS USED TO CONSERVE CORE REQUIREMENTS. 

MATRICES A, 3, C AND D CONTAIN THE RESULTS 
OF THE GLOBAL INTEGRATION PERFORMED ON EACH 
TERM IN THE EQUATION. MATRICES A AND B ARE USED 
DURING THE GAUSS-SIEDAL ITERATION ON THE U AND 
V EQUATIONS. MATRICES C AND D ARE USED DURING 
THE ITERATION ON THE GEOPOTENTIAL FIELD. 

VECTORS 

AREA2= TWICE THE AREA OF EACH TRIANGLE 
A1 =A CONSTANT OVER EACH TRIANGLE USED TO 
EVALUATE A LATITUDE DERIVATIVE. 

A2= A CONSTANT SIMILIAR IN NATURE TO A1 
A^= SAME AS Ai 

Bl= A CONSTANT USED TO EVALUATE A LONGITUDE 
DERIVATIVE. 

B2= SAME AS B1 
B3= SAME AS BI 

XLAT = LATITUDES OF ALL THE POINTS 
XLON= LONGITUDES OF ALL THE POINTS 
XLONl= LONGITUDES OF THE POINTS THAT ARE 
USED TO HAVc CYCLIC CONTINUITY 
COSLAT= COSINE OF ALL THE LATITUDES OF EACH 
POINT 

SINLAT= SINE OF ALL THE LATITUDES OF EACH 
POINT 

NTRI= NUMBER OF THE TRIANGLE WHERE CYCLIC 
CONTINUITY IS MAINTAINED THRU USE OF 
VECTOR XLCNl VICE XLON. ON ONE SIDE 
OF THE CYCLIC CONTINUITY LINE THE 
TRIANGLES HAVE 360 DEGREES AS A LONG- 
ITUDE WHILE FROM THE OTHER SIDE THE 
TRIANGLE EXPERIENCES 0 DEGREES. 

NTRO= NUMBER OF THE POINT WHERE CYCLIC 
CONTINUITY IS MAINTAINED 
TERM2= TEMPORARY HOLDING VECTOR USED DURING 
ITERATION TO MINIMIZE THE COMPUTATION 
OF THE RIGHT HAND SIDE OF THE EQUATION 
SPH= A DUMMY VECTOR USED WHEN A TRIG FUNC- 
TION IS NOT BEING INTEGRATED INTO THE 
GALERKIN SPACE BY A SUBROUTINE 
EN= THE CHANGE IN THE VARIABLE DURING ONE 

TIME STEP. THIS IS THE VECTOR THAT THE 
ITERATION SCHEME CONVERGES TO. 

NTRI5= THE NUMBER OF THE TRIANGLES THAT ARE 
SUPPORTED BY ONLY FIVE TRIANGLES VICE 
SIX TRIANGLES. THIS OCCURS AT THE 
VERTEXES OF THE MAJOR SPHERICAL TRI- 
ANGLES. 

VECT= TEMPORARY HOLDING VECTOR USED DURING 
THE PHI EQUATION. 

Ul|Vl,PHII= INITIAL FIELDS AND FIELDS AT 
THE TIME LEVEL N 

U2tV2tPHIZ= FIELDS AT TIME LEVEL N+1 
USTAR=TIME EXTRAPOLATED ASSUMPTION ON U 
• VSTAR=TIME EXTRAPOLATED ASSUMPTION ON V 
PHSTAR=TIME EXTRAPOLATED ASSUMPTION ON PHI 
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E,F AND G= HOLDING VECTORS FOR THE TANGENT, C 
CORIOLIS AND PRESSURE GRADIENT TERNS C 
PHI3= HOLDING VECTOR USED FOR HARMONIC C 

ANALYSIS C 

COA = COEFF IC lENT OF A FOR EACH TRIANGLE IN C 
EXPRESSION F=A+B*LAMBDAfC»THETA C 

COB=COEFF IC lENT FOR B FOR EACH TRIANGLE C 

BFLD=GL03AL 73 X 35 5 DEGREE GRID C 

FLD=NORTHERN HEMISPHERIC GRID USED FOR PLOT C 
CL=CGNTOUR LEVELS FOR PLOT C 

TABLES C 

NPTS=GLC3AL CORRESPONDENCE TABLE C 

NCOR=GLOBAL CORRELATION TABLE C 

NTABL=TABLE TO CONVERT FROM ICOSAHEDRAL C 

GRID TO 73 X 35 GRID. C 

COC = COEFF IC I ENT FOR C FOR EACH TRIANGLE C 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
INTEGER-2 NPT S , NCOR , NTRI 5 
INTEGER«2 NTA6L 

COMMON /CMl/ AREA2( 1296) ,A1 ( 12961 ,A2( 1296) , A3( 1296) , 
1B1(1296),B2(1296),33( 1296) 

COMMON /CM2/ XL AT ( 63^ ) , XLON ( 684 ) , XLON I ( 6 84 ) 

COMMON /CM3/ NPTS(1300,3) 

COMMON /CM4/ COSL A T ( 684 ) , S I NLAT ( 684 ) 

COMMON /CM5/ NCOR (634, 7) 

COMMON /CM6/ NTR I ( 46 ) , NTRO ( 23 ) 

COMMON /CM7/ TERM2(684) 

DIMENSION A (684, 7) ,3(684,7) ,C ( 684 , 7) , D ( 6 34 , 7) 
DIMENSION SPH(684) ,EN(684) ,NTRI5(12),VECT(684) 
DIMENSION U1 (634) ,U2 ( 684) , V 1( 684 ) , V2(o84) 

DIMENSION PHI 1(684) ,PHI2I684) 

DIMENSION USTAR( 684) ,VSTAR( 684) , PHST AR( 684) 

DIMENSION E(6S4) ,F(6o4) ,G(684) 

DIMENSION PHI3(36) 

DIMENSION NT ABL( 72,18) 

DIMENSION CO A (640) , COB (640) ,COC( o40) 

DIMENSION DFLD(73,35) ,FLD(63,63) 

DIMENSION CL (16) 

DATA OT2/360. /,DT/ 720. /,N/ 8/ , OMEGA/ 7 . 2 921 1 5432 E-05 / 
DATA WAVENO/ 8 . / , W AV V E L/ 1 . 6 1 604E -05/ , AM P / 5.49364307/ 
DATA IFLAG/0/, TIME/0 .0/ 

LOGICAL=i=l LTG( 3 ) /3=s=. FALSE ./ 

REALMS TITLEK12)/' HINSMAN',' D.E. ', 'INITIAL • 

1 , 'CONO ITIO* , 'NS WAVEN' , 'UMBER 8 ', 'PHASE SP', 

2'EED 10 D' ,' EGREES/0' , ' AY A = 7.E','+07 ', 

3' '/ 

REALMS TITLE2(12)/' HINSMAN',' D.E. ','FORCAST ' 

1 , 'CONDIT 10' , ' NS WAVEN' , 'UMBER 8 ', 'PHASE SP', 

2'EED 10 O' , • EGREES/D' , 'AY A=7.E','+07 ', 

3 ' ' / 

REAL=i''8 SU8TIT(12)/' 12',' 24',' 36', 

1 ' 48' , ' 60' , ' 72 ' , 

2' 84' , ' 96' , ' 108' , ' 120' , 

3' 132',' 144'/ 

N0TRI=1 280 
N0PTS=642 
J1=0 
JM=18 
Z = 8. 



READ(5,100) XLAT 
READ(5,100) XLON 
DO 3 I=1,N0TKI,5 
L=I+4 

3 READ( 5, 101 )( ( NPTS(K, 1 ) ,NPTS (K, 3) ,NPTS(K ,2) ) ,K = I , L) 
DO 4 I=1,N0PTS 

4 READ( 5, 102)(NC0R( I , J) ,J = 1,7) 
lUO FORMAT! 8F10. 6) 

101 FORMAT! 5(315)) 

102 F0PMAT(7I10) 

103 FORMAT! 10F8. 1 ) 

READ(5,104) NTRI5 

104 FORMAT! 1018) 
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READ(5,106) NTRI 
READ(5,106) NTRO 
106 FQRMATOIIO) 

DO 1 I=1,NQPTS 
00 1 J=l,7 
A( I , J)=0.0 
B( I , J) = 0.0 
C(I , J) = 0.0 

1 D(I,J) = 0.0 

DO 2 I=1,N0PTS 

2 SPH(I)=1. 

CALL AREAZ(iNOTRI , OMEGA , NOP T S ) 

CALL SURVEY(NOTR I ,NTABL,NOPTS) 

DO 325 1=1 ,72 
IL=73-I 

325 WRITE (6, 326) (NTABLl IL , J) , J= 1 , 1 8 ) 

326 F0RMAT(/1X,18I7) 

CALL SOLUT (XLATtXLON, Ui ,V1 , PHI 1 , NO PT S , W A VEN 0 , WAV VE L , 
1AMP,0,0,0,DT,TIME) 

U1(1)=0.0 
U1 (NCPTS)=0.0 
Vl( 1) =0.0 
VI (N0PTS)=0.0 
WRITE(6,130) 

130 F0RMAT</1X,' INITIAL PHI FIELD') 

WRITEI6, 105) ( PHI 1 ( I ) , I =1 ,NCPTS) 

WRITE(6,13i) 

131 F0RMAT(/1X,' INITIAL U FIELD') 
WRITE(6,105)(U1(I),I=1,N0PTS) 

WRITE(6,132) 

132 F0RMAT(/1X,' INITIAL V FIELD') 

WRITE(6, 105) ( Vl( I ) ,1=1 ,NOPTS) 

DO 50 1=1 , 16 

50 CL( I )=1 . 

CALL EVALINOTRI , PH I 1 , COA , COB , C OC ) 

C ALL DISPLibFLD, NTABL, COA,COB,COC) 

DO 36 J = 1 , 18 
L = 1 

JJ=36-J 

DO 37 1=2, 72,2 
PHI3I L) =BFLD( I , J J) 

37 L = L<-1 

36 CALL HARMAN! J1,JM,Z,PHI3,JJ) 

CALL LONGOU! BFLD ,FLD, PHIl ( 1 ) ) 

CALL C0NTUR(FL0,63,63,63,CL,-16, TITLE 1,6,6, LTG) 

DO 5 I=1,N0PTS 
USTARd )=U1( I) 

VSTARI I ) =V1( I) 

PHSTAR! I )=PHI1 ( I ) 

VECT( I ) =0.0 

5 EN(I)=0.0 
324 FORMAT! ////) 

KMAX=60 
DO 71 KL=1,6 
DO 70 KK = 1 ,60 
DO 7 I=1,N0PTS 
E! I) = 0.0 
F! I )=0.0 
7 G! I)=0.0 

CALL AMATRUNOTRI , A, COSLAT) 

DO 700 I=1,N0PTS 
DO 700 J=l,7 
700 C ! I , J) =A( I , J ) 

501 CALL BMATRIINQTRI ,USTAR,B1 ,B2,B3 ,SPH,A,3 ) 

CALL BMATRIiNCTRI , VST AR , A 1 , A2 , A 3 , COSL AT , A , B ) 

2000 CALL DMATRKNDTRI ,USTAk,Bl ,B2,B3,SPH,D) 

CALL DMATRI! NOTRI , VST AR , A 1 , A2 , A 3 , COS L AT , D ) 

DO 6 I=1,N0PTS 
DO 6 J=1 , 7 

C! I , J)=C! I ,J )-0T2*D!I , J) 

6 D! I ,J ) = D! I ,J ) *0T 

202 CALL SOLVE!C ,D, EN ,PHI 1 ,NOPTS,NTR I5,VECT , IFLAG ) 
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DO 8 I=l,NOPTS 
8 PHI2( I )=PHI1 ( I ) 4-EN( I) 

210 SUM=0.0 

DO 60 1=2,6 

60 SUM=SUM+PHI2 ( I ) 

AVPHI=SUM/5. 

DO 61 i =1 ,6 

61 PHI2(I)=AVPHI 
SUM=0.0 

DO 62 1=637,041 

62 SUM=SUM+PHI2 ( I ) 

AVPHI=SUM/5. 

DO 63 1=637,642 

63 PHI2(I )=AVPHI 

CALL TANVEC(NOTRI , US T AR , S I N L AT , V STAR , G ) 

204 CALL CORIO(NCTRI ,VSTAR,E,SINLAT, COSLAT) 

CALL PGF( NOTRI , PHSTAR,B1 ,B2 ,B3, SPH,F) 

DO 11 1=1 ,NOPTS 
E ( I ) = DT*( E( I)-F( I ) fG( I ) ) 

11 EN( I)=0.0 

CALL SGLVE(A,S,EN,U1,N0PTS,NTRI5,E,IFLAG) 

DO 12 I=1,NGPTS 

12 U2( I)=U1 ( I )<-EN( I ) 

201 00 13 1=1 tNDPTS 

EN( I )=0.0 
E( I)=0.0 
G( I)=0.0 

13 F ( I ) = 0. 0 

207 CALL TAUVEC(NDTRI , JSTAR,SINLAT,USTAR,G) 

CALL CORIO(NOTRI , UST AR , E , 5 I NLAT , C OSL AT ) 

CALL PGF (NOTRI , PHST AR , Al , A2 , A3 , COSLAT , F ) 

DO 15 1 = 1,N-)PTS 

15 G( I)=-DT*(E( I )+F ( I )>G( I ) ) 

2001 CALL S0LVE(A,B,EN,V1 ,N0PTS,NTRI5,G,IFLAG) 

DO 18 I=1,M0PTS 
13 V2( I )=V1( I ) +EN( I ) 

200 DO 17 I=1,N0PTS 
DO 17 J=1 , 7 
A( I, J)=0,0 
B( I , J )=0.0 
C( I , J )=0.0 
17 D( I, J)=0.0 

DO 16 I=1,N0PTS 

PHSTAR( I )=1 ,5=^PHI2 ( I ) - .5*PHI 1 ( I ) 

USTAR( I ) = 1 .5=^0 2 ( I )-.5=^Ul( I ) 

VSTAR( I ) =1.5*V2 ( i )-.5*Vl( I ) 

PHIim=PHI2(I) 

UK n=U2 ( I ) 

Vi( I }=V2( I ) 

16 EN(I)=0.0 

105 FORMAT! /IX ,6E20. 6) 

70 CONTINUE 

TIME=FLUAT( ( KL- 1 ) *KM AX +KK - 1 ) *0 T 
TIMEl=TIME/3600. 

WRITE(6,310) TIMEl 

310 FORMAT ( /IX, ’ PROG SOLUTION AT ', 2X , F5 . 1 , 2X, * HOURS ' ) 

WRITE( 6,122) 

122 F0RMAT(/1X,' PHI FIELD') 

WRITE(6 , 105) ( PHI 2 ( I ) , 1 = 1 ,NCPTS ) 

WRITE(6,110) 

no F0RMAT(/1X,' U FIELD') 

RRITE( 6, 105) ( U2( I ) ,I = 1,N0PTS) 

WRITE(6,111) 

111 F0RMAT(/1X,' V FIELD' ) 

WRITE(6,105) ( V2( I ) ,1 =l ,NQOTS) 

CALL EVAL(N0TPI,PHI1 , COA , COB , C OC ) 

CALL DISPL(3FLD,NTABL, C0A,C08,C0C) 

DO 38 J=l,18 
L=1 

JJ=56-J 

DO 39 1=2 ,72 ,2 
PHI3( L) =8FLD( I,J J ) 
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39 L=L+1 

33 CALL HAPMAN( Jl, JM,Z,PHI3, JJ) 

TITLE2( l.i)=SUbri T(KL) 

CALL LONGOU( BFLD ,FLD,PHI1 ( 1 ) ) 

CALL C0NTUR(FLU,63 ,63,63,CL,-16 , T I TL E 2 , 6 , 6 , LT G ) 
71 CONTINUE 

1^1 FCRMAT(/1X,' 3 MATRIX') 

140 FORMAT! /IX,' A MATRIX') 

107 FORMAT! /1X,7E15.6 ) 

STOP 

END 



SUBROUTINE SURVE Y ! NOTR I , NT ABL , NOPTS ) 

INTEGER^2 NPTS,NTABL 

COMMON /CM2/ XL AT ! 6 34 ) , XLON ! 684 ) , XLON L ! 6 84 ) 

COMMON /CM3/ NPFS! 1300,3) 

COMMON /CM4/ COSL A T ( 68^ ) , S I NL A T ! 634 ) 

COMMON /CM6/ NTR I !46) , NTRO! 23 ) 

DIMENSION S!3) ,ISAVE!3) ,NTA3L!72,1) 

DATA RAD/ 6. 28 31 8 53 08/, THETA /I. 57 079632 7/ 
DRAD=RA0/72. 

NOTR=NJTRI /2 
DTHET=THETA-ORAD 
NGPT=341 
DO 1 1=1, 72 
XL0=FL0AT!I-1)*0RA0 
YLA=DTHET-DRAD 
DO 2 J=2, 13 
I FLAG=0 

IF(J.EQ.18) YLA = YLA^-.01 
ISAVE! i )=1 
ISAVE! 2)=1 
40 IK0UNT=0 
L = 2 

DO 3 K=6,N0TR 

IFIIFLAG.EQ. l.AN0.K.EQ.ISAVE!2) ) GO TO 3 

IC0UNT=0 

XL0MX=0.0 

XL0MN=2.*RA0 

YLAMX=0.0 

YLAMN=2 .^RAO 

LL = 0 

IFINTRI !L ) .NE.K) GO TO 4 

L = L<-1 

LL=1 

DO 5 KL=1 ,3 

XL0MX=AMAX1!XL0MX,XLGN1 !NPTS!K,KL) ) ) 

5 XLQMN=AMI N1 !XL0MN,XL0N1 !NPTS!K,KL) ) ) 

GO TO 8 

4 DO 6 KL=1 , 3 

XL0MX=AMAX1!XL0MX,XL0N! NPTS!K,KL ) ) ) 

6 XL0MN=AMIN1 ! XLOMN , XLON ! NPTS ! K, KL ) ) ) 

8 DO 7 KL=1 , 3 

YLAMX=AMAX1 ( YLAMX , XLAT! NPTS ! K ,KL ) ) ) 

7 YLAHN=AMIM1 1 YLAMN , XLA T! NPTS 1 K , KL ) ) ) 
IF!XLO.LT.XLOMX. AND . X LO . GE . XLOMN ) ICOUN T= I COUNT^- 1 
IF ! YLA . LT .YLAMX. AND. Y LA. GE . YLAMN ) I COUNT = ICCUNT^- I 
IF! IC0UNT.lt. 2) GO TO 3 

IKCUNT= IKOUNT+1 
IF! IK0UNT.EQ.2) GO TO 9 
ISAVE! 1 )=K 
GO TO 3 

9 ISAVE!2)=K 
GO TO 10 

3 CONTINUE 
10 IC0UNT=0 

DO 11 KL=1,3 
DO 12 KK=1,3 
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IF(NPTS( ISAVE( 1) ,KL) .NE.NPTS( I SAVE(2) ,KK) ) GOTO 12 
IF( ICOUNT.GT.G ) GO TO 13 
IS1=NPTS( ISAVE( 1 ) ,KL) 

ICOUNT= ICOU^4T + l 
GO TO 11 

13 IS2=NPTS( I SAVE( 1 ) ,KL) 

GO TO 14 

12 CONTINUE 
11 CONTINUE 

14 IF(LL.EG.O) GO TO 15 
IFLAG=1 

I FIXLONU IS2 )-XLONl( I SI ) .EQ.0.0) GO TO 40 

SLOPE = ( XLAT( IS2)-XLAT( I S 1 ) J / ( XL ON 1 ( I S 2 ) - XLON 1 ( I S 1 ) ) 

GO TO 15 

15 IFLAG=1 

IF(XLON(IS2J-XLONnsi).EQ.O.O) GO TO 40 

SLOPE = lXLAT( IS2 5-XLAT ( ISl ) )/(XLON( I S 2 ) -X LON ( I S 1 ) ) 

16 IF(LL.EC.O) GO TO 17 

B = XLAT( ISn-SLOPE*XLONl( ISl) 

GO TO 18 

17 6=XLAT( ISl )-SLCPE*XLON( ISl) 

18 YLAC = SLGPE’!=XL0<-3 
DO 20 KL=1t3 

IF( ISUEQ.NPTS(ISAVE( 1) ,KU .OR.IS2.EQ.NPTS( I SAVE! 1) , 
IKD) GO TO 20 
IN0T1=NPTS( I SAVE! 1) ,KL) 

GO TO 22 
20 CONTINUE 

22 DO 23 XL=1,3 

IF( ISl.EQ.NPTSI I SAVE( 2) ,KL1.0R.IS2.EQ.NPTS(ISAVE(2), 
IKL) ) GO TO 23 
IN0T2 = NPTS(ISAVE( 2) ,KL) 

GO TO 24 

23 CONTINUE 

24 IFdNOTi .GT.IN0T2) GO TO 25 
MP=ISAVE ( 1 ) 

NM=ISAVE(2) 

GO TO 26 

25 MP=ISAV£(2) 

MM=ISAVE ( 1 ) 

26 IF(YLA.GT.YLAC) NTABL ( I , J ) =MP 
IF( YLA.LE.YLAC) NTABL<I,J)=MM 

2 YLA=YLA-DRAD 
1 CONTINUE 
DO 30 1=1,15 

30 NTABLI I , 1 )=5 
DO 31 1=16,29 

31 NTAbL ( I , 1 ) =4 
DO 32 1=30,44 

32 NTABLI I , 1 )=3 
DO 3b 1=45,58 

33 NTABLI I , 1 ) =2 
DO 34 1=59,72 

34 NTABLII ,1)=1 
NTABLI 72, 16)=559 
NTABLI 71,17)=639 
NTABLI 72, 17)=639 
NTABLI 71 ,1 8)=640 
NTABLI 72, 18)=639 
RETURN 

END 



SUBROUTINE SOLVE I A , B , EN , Z , NOPTS , NTR I 5 , C , I FLAG) 
INTEGER*2 NC0R,NTRI5 
COMMON /CM5/ NC0RI684,7) 

COMMON /CM7/ TERM2Io84) 

COMMON /CM8/ VI SI 684) 

DIMENSION AI 684, 1 ) ,BI 684, 1),ENI1),ZI1), NTR 15(1) 
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DIMENSION C( 1 ) 

DATA EPSI/.1E-05/ 

KL = 2 

NOPT=NOPTS-l 
22 LX=2 

DO 10 I=KLtNOPT 
LY = 7 

IF( I .NE.NTRISILX) ) GO TO 11 

LY=6 

LX=LX+1 

11 TERM2(I)=0.0 
DO 12 J=1,LY 

12 T£RM2( I ) = TERM2( I ) <-Z(NCCR( I,J))«B(I,J) 
10 CONTINUE 

DO 1 L=1,1C0 
ZMAX=0. 0 
ERROR=0.0 
LLL=2 

DO 2 I=KL,NOPT 
LL=7 

IF( I.NE.NTRI5(LLL) ) GO TO A 
LL = 6 

LLL=LLL+1 

4 RHS=0.0 

DO 3 J=I,LL 

IF(NCOR( I , J) .EQ. I ) GO TO 5 
TERMI=cN(NCOR(I, J ) ) *A { I , J ) 
RHS=RHS-TERM1 
GO TO 3 

5 K=J 

3 CONTINUE 

RHS=RHS+TERiM2( I ) +C ( I ) 

ZDUM=EN(NCOR( I,K>) 

ZMAX=AMAX1 (ZMAX, A8S( EN( NCOR { I , K) ) ) ) 

20 EiMlNCORI I ,l<n=RMS/A( I ,K) 
IFUFLAG.NE.Z) GO TO 21 
EN(NCOR{ I ,K) ) = EN( NCDR( I , K ) ) ’i'V I S ( I ) 

21 DELTA=ABS( ZDUM-E N ( NCOR ( I ,K) ) ) 

2 ERRGR=AMAXl( ERROR, DELTA) 

EPS=ZMAX^EPSI 
IF(ERROR.LE.EPS) GO TO 6 
1 CONTINUE 

6 WRITE(6,100) L 

100 FORMAT (/IX, I 4, 2X, ’ ITERATIONS' ) 
IFIL.EQ.IOO) IFLAG=3 
RETURN 
END 



SUBROUTINE BMATRI (NOTRI ,DUM,C1 , C 2 ,C3 , S PH , A , B ) 

INTEGER=i'2 NRTS 

COMMON /CM3/ NPTS( 1300,3) 

DIMENSION DUM(1 ) ,C1(1),C2{1),C3(1),SPH(1),A(684,1), 
16(684,1 ) 

DATA EARTH/6.371 E+06/ ,DT2/360. /, DT/ 720./ 

DO 1 1=1, NOTRI 
II=NPTS( I ,1) 

JJ = NPTS( I ,2) 

KK=NPTS (1,3) 

CALL SEARCH(II,JJ,KK,I1,I2,13,II) 

CALL SEARCHl iI,JJ,KK,Jl,J2,J3,JJ) 

CALL SEARCH( I I , J J,KK,K1 ,K2, K3,KK) 

FAC1=1./(EARTH*120.) 

TERMA = 6. *DUM( I I ) ^’'^ S PH ( I I ) *■ 2 . *DUM( I I )*SPH( JJ ) 

TERMB = 2 .«DUM( I I ) =^SPH( KK ) f 2 . «OJM ( J J ) *SPH ( I I ) 

TERMC = 2 .’!=DUM( JJ ) - S PH ( J J ) +■ DUM ( J J ) * S PH ( RK ) 

TERM0 = 2.=^DUM(KK) ^SPH( I I ) +OUM ( KK ) *S PH ( J J ) 1 2 . *DUM ( KK ) * SP 
1 H ( KK ) 

TERM1= ( TERMA+TERMB+TERMC+TERMD) *FACT 
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T ERM6 = 2 .*DUM (II)«SPH(II)+-2. *OJM( 1 1) *SP H ( J J ) *-DUM ( I I) SP 
1H( KK) 

TERrtC=2 .*DUM( JJ) *SPH( I I )4-6. *DUM( JJ )*SPH( J J ) +-2.*DUM( J J) 
1«SPH(KK) 

TtRMD=OUM( KK )«SPH( I I ) 4-2 . *DUM ( KK ) « SPH ( JJ ) +2 . '-J'DUM ( KK ) * SP 
IH (KK) 

TERM2=( TERH3+TERMC +-T ER MD ) *F ACT 

T£RMC = 2 .*DUM ( 1 1 ) #SPri { I I ) +-01JM( I I ) *SPH( JJ ) +-2. «OUM( I I ) 
IH(KK) 

TERMD=DUM( J J ) «S PH ( I I ) +2 .*DU M ( J J ) « S PH ( J J ) +2 .*OUM( J J )*SP 
IH(KK) 

TERME = 2 .«DUM(KK) *SPH( I I ) + 2 . *DUM ( KK ) T- SPH ( J J ) +6 ( KK) 

1«SPH(KK) 

TERM3=( TERMC+-TERMD4-TERME)T-FACT 
A(II,Ii)=A(II,Il ) +-0T2*TERMi=:‘Cl { I ) 

6(II,I1)=B(II,I1 )-OT«TERMl^.'-Cl( I ) 

A( 1 1 , 1 2 ) = A( I I , 12 ) *-DT2'-:=TERMl *C2( I ) 

B( II ,12 )=B( 1 1 , 12 )-0T=i'TEKMl*C2( I ) 

A( II , 13 ) = A( I I , 13 )+0T2*TERMl*L3( I ) 

B(1I,I3)=3(I1,I3 )-DT*TEkMr-.--C3( I ) 

A(JJ,J1)=A(JJ,J1) +DT2>!'T£RM2’?C1 ( I ) 

B(JJ,Ji)=B(JJ,Jl )-DT*TERM2^Cl ( I ) 

A( JJ,J2) =A(JJ,J2)+DT2^TERH2*C2( I ) 

B( JJ, J2)=B( JJ, J2 )-DT=:=TERM2*C2( I ) 

A( JJ , J3 ) =A( JJ , J3 ) +0T2*TERM2*C3 ( I ) 

B( JJ, J3 ) =B( J J , J3 ) -DT-:=TERM2^-C3 ( I ) 

A(KK,K1) = A(KK,K1 ) +-DT2 *TERM3=i=Ci ( I ) 

6(KK,K1 )=3(KK,K1 ) - OT«T ERM 3*C 1 ( I ) 
A(KK,K2)=A(KK,K2)+DT2*TERM3*C2( I ) 

B(KK,K2) = 3(KK,K2 )-OT'TTERM3=^C2 (I ) 

A(KK,K3)=A(KK,K3 ) 4-DT2*TERH2^C3 ( I ) 

1 8(KK,K3)=B(KK,K3 )-DT=i'TERM3^C3 ( I ) 

RETURN 

END 



SUBROUTINE DMATRKNOTRI , DUM , C 1 , C 2 ,C3 , S PH , D ) 

INTEGER=!=2 NPTS 

CCMMGN /CM3/ NPTS (1300,3 ) 

DIMENSION DJM (1) , Cl( 1) ,C2( 1 ) ,C3( 1) ,SPH( 1 ) ,D (684, 1 ) 
DATA EARTH/o.371E^-06/ 

DO 1 I=1,N0TRI 
I I=NPTS (1,1) 

J J=NPTS (1,2) 

KK=NPTS( I ,3) 

CALL SEARCH( 1 1 , J J , KK , 1 1 , I 2 , I 3 , I I ) 

CALL SEARCH( I I , J J , KK , J 1 , J2 , J3 , J J ) 

CALL SEARCH( I I , J J , KK , K 1 , K2 , K3 , KK ) 
FACT=1./(EARTH*120.) 

TA=SPH( I I )^DUM( I 1 ) 

TB=SPH( JJ)*DUM( I I ) 

TC = SPH( KK) =^-DUM( I I ) 

TD=SPH( I I )*DUM( J J ) 

TE=SPH( JJ)«DUM( J J) 

TF = SPH( KK)*DUM( J J) 

TG=SPH( I I )^DUM(KK) 

TH=SPH( JJ)''^-DUM(KK) 

TI =SPH( KK )^=DUM( KK ) 

TERM1=6 .*TA + 2 .-( rB+-TC + TD+rE + TG«-TI) + TF + TH 
TERM1=TERMK-FACT 

TERM2 = 6 .-TE + 2,-T( TA + TB+TD+-TF +TH4-TD+TC+-TG 
TERM2=TERM2*FACT 

TERM3 = 6.#TI <-2.-( TA +-TC +• TE+- TF +TG+-TH ) +TB + TD 

TERM3 = TERM3=!^FACT 

D( 1 1 , II ) = D( 1 1 , 1 1 ) +-TERM1*C1 ( I ) 

D( I I , 12 ) =D( I I , I 2 )<-TERM2*Cl ( I ) 

D( I I , 13 )=0( 1 I, 13 ) <-TERM3*CU I ) 

D(JJ,J1)=D( JJ,J1 )+TERMl^C2( I) 

D( JJ , J2 )=D( JJ , J2 ) 4-TERM2=:-C2 ( I ) 
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D( JJ, J3)=0(JJ, J5)+TERM3*C2( I) 
0(KK,K1) = 0(KK,K1 ) + TE RM 1 *C3 ( I ) 
D(KK ,K2 )=D( KK, K2 ) +TERM2*C3 ( I ) 
1 D(KK,K3 )=D(KK,K3 )+TERM3’i=C3( n 
RETURN 
END 



SUBROUTINE AREAZ ( NOTRI , OMEG A , NOP TS ) 

INTEGER=^'2 NPTS 

COMMON /CMl/ AREA2(1296) , A1 ( 129 6) t A2( 1296) , A3( 1296) , 
1B1(L296),32( 1296) ,B3( 1296) 

COMMON /CM2/ XLAT(634) ,XLON(684) ,XL0N1(684) 

COMMON /CMji/ NPTS ( 1300,3) 

COMMON /CM4/ COS L AT ( o 84 ) , S I NL AT ( 684 ) 

COMMON /CM6/ NT R I ( 46 ) , NTRO ( 23 ) 

DATA RAD/6. 283185308/, L/1/, LL/1/ 

DO 6 I=1,NQPTS 
6 XLON( I ) =RAO-XLON( I ) 

DO 2 J=1,N0PTS 
IF(L.GT.23) GO TO 3 
IF(NTRO(U .NE. J) GO TO 3 
L = L+1 

XLONl ( J)=XLON( J) -RAD 
GO TO 8 

3 XLONK J)=XLON( J) 

8 COSLAT( J)=COS(XLAT( J) ) 

SINLAT( J)=SIN(XLAT( J) ) 

2 CONTINUE 
XLGNI ( 1 )=0.0 
XLONl(NOPTS)=0.O 
DO 1 1=1 , NOTRI 
I I =NPTS (1,1) 

J J=NPTS( I ,2) 

KK=NPTS (1,3) 

IF(LL .GT.46) GO TO 4 
IF(NTRI(LL).NE.I ) GO TO 4 
LL=LL+1 

Al( I )=XL0N1(KK)-XLCNI( JJ) 

A2( I )=XL0N1( I I )-XLONl( KK) 

A3(I)=XL0N1( JJ)-XLCN1( II) 

GO TO 5 

4 Al( I ) =XLON(KK) -XLON( J J ) 

A2( I ) =XLON( I I )-XLON(KK ) 

A3( I)=XLON(JJ)-XLON( II ) 

5 Bid )=XLAT( JJ)-XLAT(KK) 

B2( I ) =XLAT(KK)-XLAT( I I ) 

1 B3( I)=XLAT( I I)-XLAT(JJ) 

N0TR=N0TRI-4 
DO 10 1=1,5 
A2( I )=0.0 

10 A3(I)=0.0 

DO 11 I=NOTR, NOTRI 
Al( I )=0.0 

11 A2( I )=0.0 

DC 12 I =1 , NOTRI 

12 AREA2(I ) = ABS(A1( I )*62(I)-A2( I)=.'=B1( D) 

DO 14 I=NGTR,NOTRI 

14 AREA2 ( I )=AREA2( 1 ) 

DC 15 1=1,5 

15 A3( I) =-Al( I) 

DO 16 I=NUTR, NOTRI 

16 A1(I)=-A3(I) 

RETURN 

END 
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SLBROUT INE S E A^C H ( I , J , K , 1 NDEX 1 , I NDEX2 , I NDEX3, I I ) 

INTEGER*? NCGR 

COMMON /CMb/ NC0R(684»7) 

DIMENSION INDEX! 3) 

JJ = I 

DO I M=l,3 
IF(M.EQ.3) JJ=K 
DO 2 L=l,7 

IF(NCOR( II ,L I.NE.JJ) GO TO 2 
INDEX! M )=L 
JJ=J 
GO TO I 
2 CONTINUE 
1 CONTINUE 

IND6X1=INDEX!I) 

INDEX2= INDEX !2 ) 

IN0EX3= INDEX!3) 

RETURN 

END 



SUBROUTINE AMATR I ( NOT RI , A , DUM ) 

INTEGER*2 NPTS 

COMMON /CMI/ AREA2(1296) ,A1 ! 1296) ,A2! 1296) , A3! 1296) , 
IBl ! 1296) ,B2! 1296) ,B3! 1296) 

CCMMCN /CM3/ NPTS!1300,3) 

DIMENSION A! 684, 1 ) ,OUM! 1 ) 

DATA EARTH/o.371 E + 06/ 

DC 1 I=l,NOTRI 
II=NPTS! 1,1) 

J J=NPTS ! I ,2) 

KK=NPTS! I ,3) 

CALL S£ARCH!II, JJ,KK, 11,12,13,11) 

CALL SEARCH! I I , J J , KK , J 1 , J2 , J3 , J J ) 

CALL SEARCH! I I , J J , KK , K 1 , K2 , K3 , KK ) 

FACT=AREA2!I )/120. 

TERM1=6.*0UM! 1 1 ) +2 .*OUM! J J ) <-2.*DUM!KK ) 

TERH2 = 2,*DUM! ID i-2 .*DUM! J J ) +DJM! KK) 

TERM3 = 2 .*DUr-U 1 1 ) ^■DUM1 J J ) *-2 .*DUM! KK) 

TERM4=2 ,*OUM! 1 1 ) +6 .*DUM! JJ ) +2. *OUM!KK) 

TERM5 = DUM! II )^-2 .*OUM! J J )+-2. *DUM! KK ) 

TERM6 = 2 .*DUM! II ) +2 .*OUM !J J ) +6. *DUM IKK ) 

A{II,I1) = A!II,I1) 4-FACT*TERM1 
A!II,I2) = A!II,I2) <-FACT*TERM2 
A! II, 13) = A! I I, I5)+FACT*TERM3 
A(JJ, J1)=A! JJ, Ji ) +FACT*TERM2 
A! JJ, J2)=A! JJ, J2 )+FACT*TERM4 
A! JJ, J3)=A!JJ, J3 )+FACT*TERM5 
A! KK,K1 ) = A!KK,K1 ) <-FACT*TERM3 
A!KK,K2)=A!KK,K2 )+FACT*TERM5 
1 A!KK,K3)=A!KK,K3 )+FACT*TERM6 
RETURN 
END 



SUBROUTINE CORIOINOTRI , DUM, VEC T , SPHl , SPH2 ) 

IMPLICIT REAL*8 !T) 

INT£GER*2 NPTS 

COMMON /CMI/ AREA2 !1296) ,A1 ! 1296) ,A2! 1296) , A3I1296) , 
IBl ! 1296) , B2! 1296) , 33! 1296) 

COMMON /CM3/ NPTS! 1300,3) 

DIMENSION DUM!1 ) , VECT ! 1 ) , SPHl ! 1 ) , SPH2! 1 ) 

DATA EARTH/6.3 71 E + Ob/ , GMEGA/7 . 29 21 15432E-05/ 

DO 1 I=1,N0TRI 
I I=NPTS! 1,1) 

JJ=NPTS ! I ,2) 
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KK=NPTS ( 1,3) 

FACT=AREA2( I )«OMEGA/360. 

TAAA=DUM( II)^SPHl(II) =:=SPH2 (II) 

T AAB=DUM( 1 1 )’i‘SPHi ( I I )’!'SPh2( JJ ) 

TAAC=DUM( II )*SPH1 ( II )=!=SPH2( KK) 

T ABA=DUM( in^S PH 1(JJ)=:^SPH2( II) 

TABB=OUM( II ) -SPHl ( J J ) =!-- S PH2 ( J J ) 

TABC = DUM( 1 1 ) =:=SPH 1 ( J J ) *SPH2( KK) 

TACA=DUM( 1 1) «S PH 1(KK)*S?H2( II) 

TACB=DUM( II) =4=5 PHI ( KK ) =4= S PH2 ( J J ) 

TACC = OUM( II ) =4=SPH 1 ( KK) =4=SPH2 { KK ) 

TBAA = DUM( JJ) =4=SPHi ( I I ) *SPH2( I I ) 

TBAB = OLIM( J J) =.= SPH1 ( II )«SPH2( JJ) 

TBAC = DUM( J J )--!=SPHl ( I I )=4=SPH2 ( KK) 

TBBA = OUM( JJ)-SPH1 ( JJ)=v=SPH2( I I ) 

TBBB=DUM( J J) *SPH1 ( JJ ) *SPH2 ( J J ) 

TBBC = DUM( JJ )vSPHl ( J J ) S PH 2 ( KK ) 

TBCA = DUiM( JJ)=!=SPH1 (KK)*SPH2{ II) 

T6CB=DUM( JJ)^SPH1 (KK)-SPH2 ( JJ) 

TBCC = DUM( J J) =i=SPHl ( KK) *SPH2( KK) 

TCAA=DUM(KKJ*SPH1 ( II )=4=SPH2( I I ) 

TCA6 = DUM( KK)=!'-SPH1 ( I I )=!=SPH2 ( JJ ) 

TCAC=OUM( KK) =i=SPHl ( II )*SPH2(KK) 

TCBA = DUM(KK)=:-'SPH1 ( J J ) =4=SPH2 ( I I ) 

TCBB=DUM( KK)*SPH1 ( JJ)*SpH2( JJ) 

TCBC=DUM(KK)-SPH1 ( J J ) ='^SPH2 ( KK ) 

TCCA = DUM( KK)=4=SPH1 (KK)=4^SPH2( II) 

TCCB=DUM( KK)=4=SPHi ( KK ) =4=SPH2 ( J J ) 

TCCC = DUM(KK)*SPH1 ( KK ) *S PH 2 ( KK ) 

TERMA=6 .*( TAA6+T AAC + T ABA + TACA+T3AA + TB3B^-TCAA^■TCCC ) 
TERM6 = 4 .*( TABBHACC + TBAB+TB3A<-TCAC+TCCA ) 

TERMC=2 .=;= ( TABC+-TACB4-TBAC + TBBC + TBCA + TBC8+-TBCC) 
T£RHC=TERMC + 2.=4=( TCAB + TC6A.4-TC66+-TCBC«-TCCB ) 

TERMK = 24.*TAAA + TERMA+TERMB + T ERMC 
TERMK=TERMK^=FACT 

TERMO = 6.=4= ( TAAA + TABB*•TBAB+T6BA^■T8bC^-TBC3^-TCBB^■TCCC) 
TERHc=4.»M TAAB + TABA + TBAA^-TSCC + TCBC+TCC0 ) 

TERMF = 2 .=;=( TAAC+TABC+TACA + TACB + TACC<-TBAC i-TBCA) 

TERMF = T£RMF<-2.=4=( TCAA + T CAi3f TCAC+ TCBA vTCCA ) 

TERMH=24.*TBB3+TERMCtTERME+T£RMF 

TERMM=TERMM-FACT 

TERMG = 6 (TAAA + TACC + TB3B + TBCC + TCAC+-TCPC4-TCCA+TCCB ) 
TERMH=4 .=^'{ TAAC + TAC A + TB8C^■TBC6^-TCAA+-TCB8) 

TERMI = 2 .*{ TAAB^■TAbA + TA3B + TA3C^■TACB^-TBAA+TBAB+TBAC) 
TERKI = TERMl+-2 .=!=( TBBA + TBCA<-TCAB+TC3A) 

TERM0 = 24.--4=TCCC + 7 ERMG-t-TERMH + TERM I 
TERMO=TERHO--:=FACT 
VECT( I ! ) =VECT( 1 1 ) + TERMK 
V£CT( J J )=VECT( JJ )+TERMM 
1 VECT(KK)=V£CT(KK)+TERMO 
RETURN 
END 



SUBROUTINE PGF(NGTRI , DUM , C 1 , C2 , C 3 , S PH, F ) 
IMPLICIT REAL4=8 (T) 

INTEGER=4'2 NPTS 

COMMON /CM3/ NPTS( 1300,3) 

DIME NS I ON OJM(i ) , Cl ( 1 ) , C2 ( 1 ) , C3 ( 1 ) , S PH ( i ) , F ( 1 ) 
DATA EARTH/6. 371E^-06/ 

DO 1 1=1,NGTRI 
II=NPTS( 1,1) 

JJ=NPTS (1,2) 

KK=NPTS ( I ,3) 

FACT=1 ./ ( EARTH*24. ) 

TERMA = DUM( 1 1 ) *C I ( I ) ===5 PH ( I I ) 

TERMB = DUM( 1 1 ) <=C 1 I I ) =i=S PH ( J J ) 

TERNC = DUM( 1 1 ) =!=C 1 ( I ) =<=SPH ( KK) 

TERMO = DUM< JJ ) =4'C 2 ( I ) =4=S PH ( I I ) 
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TEP.ME = DUM( JJ)«C2( I )*SPH( JJ) 

TcR.MF = DUM( J J ) *C 2 { I ) *S PH ( KK ) 

TERMG = DUM( KK) *C3 ( I ) PH { I I ) 

TERMH = DUM( KK)>)'C3 ( I )«5PH< J J) 

TERMI=UUM(KK)«C3( I )*SPH(KK) 

TERMJ = 2 TERMA + TERMD +TERMG ) 

TERMK = TtRMJ+TERMB*-TERMC+TERME+TERMF<-TcRMH + TERMI 

TERMK = TEPMK>:=FACT 

TERML=2.*( TFRMB+TERME +TERMH) 

TERh^^ = TERML4•TERVlA^■TERMC^-TERMDfTERMF + TERMG^-TERMI 

TEPMH=TERMM*FACT 

TERMN = 2 .*( TERMC<-TERMF +TERMI ) 

TERMO=TERMN + TERMA<-TERMB + TERMD + TERME*-TERMG + TERMH 
TERMO=TEkMO*FACT 
F ( II )=F( I I ) tTERMK 
FUJJ=F( JJ )fTEPi^M 
1 F(KK)=F(KK)<-TERMO 
RETURN 
END 



SUBROUTINE T ANVE C ( NOT R I , OU M , SPH 1 ,SPH2,G) 

IMPLICIT REAL*8 (T) 

INTEGER^Z NPTS 

COMMON /CMl/ AREA2(1296) , A1 ( 1296) ,A2( 1296) , A3( 1296) , 
181(1296), 3 2( 129 6), B3 (1296) 

COMMON /CM3/ NPTS (1300,3) 

DIMENSION DUH( 1 ) ,SPH1 ( 1) ,SPH2( 1 ) ,G(1 ) 

DATA EARTH/6. 37ie<-06/ 

DO 1 I=1,N0TRI 
II=NPTS(I, 1) 

J J=NPTS (1,2) 

KK=NPTS (1,3) 

FACT=AREA2( I )/( 7 20 . ’i^EA RTH ) 

TAAA = DUM( ID^^SPHK II)*SPH2( II) 

TAAB=DUM( ID^i^SPHK II)=:'SPH2( JJ) 

TAAC=DUM( II ) ^i'SPHl ( I I )^SPH2 ( KK) 

TABA=DUt«l( II )*SPH1 ( JJ )*SPH2( I I ) 

TABB=DUM( II ) *SPrii ( JJ )=i'-SPH2 ( JJ ) 

TABC = DUM( II)-SPH1 ( JJ)*SPH2(KK) 

T ACA=DUM( II ; *SPH1 ( KK) *SPH2 (II) 

TACB=DUM( II)-SPH1(KK)*SPH2( JJ) 

TACC = DUM( II)*SPH1 ( KK ) *SPH2 ( KK ) 

TBAA = DUM( J J )*SPH1 ( I I )-^SPH2( I I ) 

TBAB=DUM( JJ)*SPHI ( II )*SPH2( JJ) 

TBAC = DUM( JJ )=!^SPH1 ( I I )*SPH2( KK) 

TBBA = DUM( J J) ^.= SPH1 ( JJ )-‘'SPH2 (II) 

TBBB = DUM( JJ)*SPH1 ( JJ)*SPH2( JJ) 

TbBC=DUM( J J )=:=SPril ( J J ) *SPH2( KK) 

TBCA=DUM( JJ)'i‘SPHl ( KK)*SPH2( I I ) 

TBCB = DUM( J J ) vSPHl ( KK ) *SPH2 ( J J ) 

TBCC=DUM( JJ)^'SPH1 ( KK ) S PH2 ( KK) 

TCAA = DUM( KK)*SPH1 (II) '-=SPH2 (II) 

TCAB=DUM( KK)*SPHi ( I I )=i=SPH2( JJ) 

TCAC = OUM( KK)«SPH1 ( I I ) *S PH2 ( KK ) 

TCBA=DUM( KK)*SPH1 ( JJ )*SPH2 ( II) 

TCB6=DUM( KK)=:‘SPH1 ( JJ)*SPH2( JJ) 

TCBC = DUM( KK) -SPHi ( J J ) * SPH 2 ( KK ) 

TCCA = DUM( KK) ^i'SPHl ( KK) -SPH2 (II) 

TCCB = DUM( KK)*S?H1 (KK)«SPH2( JJ) 

TCCC=DUM( KK ) *SPH 1 ( KK) «SPH2 ( KK ) 

TERMA = 6 TA AB + TAAC + T ABA + TAC A+T8AAI-TB8B + TCAAf TCCC ) 

TERMB = 4.v( TABB + TACC<-TBAB^-TBBAfTCAC<-TCCA ) 

T£PMC = 2.*(TAbC»-TACB + TBAC + TBBC + TBCA<-TBCB<-TBCC) 

TERHC = TERMC f2.*(TCAB<-TCBA4-TCBrt + TCBC<-TCC3) 

TERMK = 2 4. *TAAA + T ERMA*- TERMB^- TER MC 
TERMK=TERMK-FACT 

TERMD = 6 ( TAAA^■TABB^■TBAO^■TeBAf T3BCfTBCa + TCBB<-TCCC ) 

TERME=4 .=!= ( TAAB<-TABA<-TBAA + TbCC ♦•TCBC + TCCB ) 
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TERMF = 2 .*( TAAC + TABC+TACAtTACB«-TACC+TBAC + TBCA) 

TERMF = TERMF + 2 TCAAfTCAB + TCAC^-TCBA^-TCCA) 

TERMM = 2 4.=i'T386+TERMJ+TERME + TERMF 
TERMM = TE RMM"=F ACT 

TERMG=6 .=!= ( TAAA^•TACC + TBB6^-T8CC^■TCAC^■TCBCf TCCA+-TCCB ) 
TERMH=4.*( TAAC+TACA+TBBC+TBCB+TCAA+TCBB ) 

TERMI = 2 ( TAAdi-TAB A4-T ABB<-TABC + TACB + TBAA + TBAB+TBAC ) 

TERNiI =TERMI+2.*( TB6A + TBCA+ TCAB + TCBA) 
TERM0=24.*TCC:+TERMG+TERMH+TERMI 
TERFiO = TERMOvFAC T 
G( II ) = G( I I J+TERMK 
G( Jj)=G( JJ J+TERMM 
1 G(KK}=G(KK)+TERMO 
RETURN 
END 



SUBROUTINE SOLUTIXLATI ,XLONl,UI , VI , PHI I , NOPTS i WAVENO , 

1 WAVVEL, AMP,KL,KK,KM,DT ,TIME) 

DIMENSION XL ATI ( I ) ,XLONl( 1 ) , U I ( 1 ) , V I ( 1 ) , PH I H I ) 

DATA EARTHA/6.37I E+06/tOMtGA/7.29 211i)43 25E-05/ 

DATA BASEHT/5570./,GRAV/9.8/ 

ASTAR=AMP«EARTHA 
WNPl =WAVENO+I , 

PH AZSP= WAVVEL/ HAVE NO 

B=(PHAZSP+2.-OMeGA/(WNPl«(WMPH-l.) =5= ( 1 . / ( 1 .-2 . / ( WNP 1* 

1 (wNPi + 1 .) n ) 

A2=EARTHA«^-'2 
DO 1 I=l,NO?TS 

TERM1 = d/ 2.'^( 2.*OMEGAf B)*(COS(XLATl (I ) }** 2 ) 

TERM2A= .25-'( ASTAR/IA2) 

TERM2B=COSl XLATI ( I ) (2* I F IX( WAVENO) ) 

IF (TERM2B.lt. .IE- 20) TERM2B=0.0 
TERM2C = WNPl"'COS( XLATI ( I ) ) *«2 
TERM2D = 2 .-''WAVENO*^:' 2- WAVENO- 2. 

TERM2E = 2.--'WAVEN0'A'*2/ (C0S(XLAT1( I ) )«’i^2) 

ATHETA = TERM1 +TE R M2 A=:=TER M2 B* ( TERM2C + TERM2 D- T ERM2 E ) 

TERM! =2 .*( OMEGA+b) -ASTAR/A2 
TERM2 = WNP1«( WNPU-1 . ) 

TERM3*=C OS ( XLATI ( I ) ) -* I F I X ( W AV ENO ) 

TERH4 = WAVEN'J*=:'2t2.^^WAV£N0t2. 

TERM5 = WNPP>‘-2*CQS( XLAT I( I ) )^^=2 
BTHETA = TERMl/TERM2-TERM3^.'-( TERM4-TERM5 ) 

TERM1 = . 25=i=(ASTAR/( A2) )-''-=:=2 

TERM2 = COS( XLATI ( I ) ) «- ( 2=^ I F i X ( WA V ENO ) ) 

IF(TERM2.LT. .IE-20) TERM2=0.0 
TERM3=WNPr-.‘'COS( XLATI ( I ) ) >^=*2- ( W N P 1 + I . ) 
CTHETA=TERMI*TERM2«TERM3 
TERM1=A2*ATHETA 

TERM2 = A2*BTHETA4^S INIXLONl ( I ) « W A V ENO- WAV V EL« T I ME ) 

TERM3 = A2=i‘CTHETA* ( 2 .*S I N( XLONl ( I ) ^WAVENO- WAV VEL-T I ME ) 
l’i^=i'2-l . ) 

PHII ( I > =TEkMH-TERM 2 + TERM3<-8ASEHT*GRAV 
TERM1=ASTAR-SIN( XLONl ( I ) *WAVENO- WAVVEL«T IME) 

TERM2 = COS(XLATl ( I ) ) I F I X ( WN PI ) 

TERM3 = WAVEN0=!^ASTAR«SI N( WA VENO=^XLONl ( I ) - WAVVEL^^'T I ME ) 
TERM4=C0S( XLATI ( I ) ) «=:=I F I X ( WAVENO- 1. ) 

TERM5 = SIN( XLATI ( I ) )**2 
TERM6 = B’i^A2-C0S( XLATI ( I ) ) 

UI ( I ) =-l./EARTHA*( TERM 1=:'TERM2-T ERM3- TER M4*TERM5-TERM6) 
TERM1 = ASTAR=^WAVEN0*SIN(XLAT 1(1)) 

TERM2 = COS(XLATl( 1 ) ) *- I F I X ( W A V ENO - 1 . ) 

TERM3 = C0S( XLONl ( I ) ^^‘WAVENO- WAVVE L^l-'T I ME ) 

I VI( I)=i ./EARTHA*( TERM 1=!' TE RM 2*T E RM3 ) 

RETURN 

END 
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SUBROUTINE EVAL( NOTRI ,DUM, A , B, C ) 

INTEGER«2 NPTS 

COMMON /CM2/ XL A T ( 6 84 ) , XLON ( 684 ) , XLON I ( 6 84 ) 
COMMON /CM3/ NPTS( 1300,3) 

COMMON /CM6/ NT R I ( 46 ) , NTRO ( 23 ) 

DIMENSION OUM( 1 ) , A( 640) ,B{ 640) ,C (640) 

N0TR=N0TRI/2 
DO 1 I=l,NOTR 
II=NPTS( I ,1) 

J J=MPTS( I ,2) 

KK=NPTS (1,3) 

TA=DUM( II)-DUM( JJ) 

IF(NTRI ( L) .NE. I ) GO TO 2 
TB=XLONl ( ID-XLJNl ( KK ) 

TD=XL0N1( ID -XLONl ( JJ) 

GO TO 3 

2 TB=XLON( I I)-XLON(KK) 

TD=XLON( I I ) -XLQN ( J J J 

3 TC=OUM( II )-DUM(KK) 

TF=XLAT( I I )-XLAT ( J J) 

TG=XLAT( I I )-XLAT(KK) 
C(I)=(TA*TB-TC^TO)/(TB«TF-TG*TD) 

IF(TD.EQ.O) GO TO 5 

B(I) = (TA-C(I )*TF)/TD 
GO TO 6 

5 B(D = (TC-C(I)=5=TGJ/TB 

6 I F(NTRI ( L) .NE. I ) GO TO 4 

A( I ) = DUM( II )-C ( I )«XLAT( II )-B( I )=l=XLONl( I I ) 

L = L+1 
GO TO 1 

4 A( I ) = DUM( ID -: ( I )*XLAT( I D -B( I )*XLON( I I ) 

1 CONTINUE 

RETURN 

END 



SUBROUTINE D I SPL ( BFLD , NTABL , COA , COB, COG ) 
INTEGER*2 NTABL 

DIMENSION BFLD (73, 35) , NT A3L ( 72 , 1 8 ) 

DIMENSION COA(l) , C03( 1) ,COC( 1) 

DATA RAD/6.283185308/, THETA/ 1. 570796327/ 
DRAD=RAD/72. 

DTHET=THETA-DRAD 

YLA=DTHET 

DO 1 J=1 , 18 

DO 2 1=1,72 

XLO=FLOATn-l)*DRAD 

TA=COA( NTABL( I , J ) ) 

TB=COB( NTABL( I , J ) ) 

TC=COC( NTABL ( I , J ) ) 

J J=36-J 

8FLD( I, JJ)=TAfTB«XLO+TC*YLA 

2 CONTINUE 

1 YLA=YLA-DRAD 
DC 3 J=18,35 

3 BFLD( 73 , J) =BFLD( 1 , J) 

TD=-. 2617993378 
TE=-. 1745329252 
TF=-. 0872664626 

BFLD( 72,2D=C0A( 430)fCOB(48O)*TF-COC(4SO)«TD 
BFLD( 72 , 20) = C0A( 559) +COB ( 55 9 ) =;DF-C0C ( 55 9) =:=TE 
BFLD( 71 , 20) = C0A( 5 60) +COB ( 560 ) ^TE-COC ( 560 ) =!^T E 
3FLD( 7 2, 19) = C0A( 63 3) +COC ( 6 3 8 ) •:= T F-COC ( 6 8 ) -T F 
BFLD( 71 , 19) = C0A( o39) 4-COB ( 63 9 ) * T E- COC ( 63 9 ) «T F 
BFLD( 7 2,18)=C0A( 63 8) 4-C0B(638)=^TF 
BFLD(71,18)=COA(64 0)4-COB(690)*TE 
BFLD( 70 ,1 8 )=C0A( o4 0) 4-COB(64G)«TO 
SUM=0.0 
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DO 5 1=1,72 
5 SUM=SUM4-BFL0( 1,13) 
AVPHI=SUM/72 . 

DO 4 J=1 , 1 7 
DO 4 1 = 1 , 73 
4 BFLC( I , J) =AVPHI 
RETURN 
END 



SUBROUTINE L0N30U ( BFLD, FLO , XNOPO) 
DIMENSION FLD(63, 63) ,BFLD( 73,35) 

DO 10 1=1,63 
DO 10 J=1 ,63 
AJ=J-1. 

AI=I-1 . 

CALL LLCVT3IAI ,AJ,ALAT, ALON) 

BI=( 3 50.0-AL0N)/5.0-H.0 
IF(BI.LT.l.O) BI=(710.0-AL ON )/5. 0 + 1.0 
I F( ALAT ,GT .84.99 ) ALAT = S4.99 
BJ=( ALAT/5.0)+18.0 

10 CALL CINTRPC BI ,BJ , BFLD, FLO( I, J) ,73,35) 
FLD( 31 ,31 )=XNOPO 
RETURN 
END 



SUBROUTINE CINTRPI FI I , F J J , F , AFF , K , L ) 

DIMENSION FIK,L) 

I = Fi I 
J = FJJ 
AI = I 
BJ = J 

R=FII-AI 
S=FJJ-B J 

IF(I.EO.l) GO TO 10 
I F( I . EQ .K-1 ) 30 TO 10 
IFI J.EQ.l ) GO TO 10 
IFU.EQ.L-l ) GO TO 10 

A=( ( F( I , J + 1 )-F( I , J ) ) f ( F(I , J ) - F( I , J-1) ) )*.5 
B=3.-(F( I,J+l)-F(I,J))-(2.*A+((F(I,J+2)-F(I,J+l))f(F(I 
1, J + l)-F( I,J) ) )*.5) 

C=( A+ ( ( F( I , J + 2)-F( I, J+1 ) ) + ( F( I , J + l)-F( I, J ) ) .=■•= (F 

1(I,J+1)-F( 1,J)) 

F1 = F( I, J ) +S*( A + S*! B+S*C) ) 

1 = 1+1 

A=( ( F( I , J+1 )-F( I , J) ) + ( F( I , J )-F ( I , J-1 ) ) )*.5 
B=3.*(F( I , J+l ) -F ( I , J ) ) - (2 .^i'A+I ( F( I , J+2) - F( I , J + 1 ) ) + ( F( I 
1, J + l)-F( I, J) ) )«.5) 

C=(A+((F(I,J+2)-F(I,J+l))+(F(I,J+l)-F(I,J)))*.5)-2.*( 
1F( I, J + l )-F( I,J) ) 

F2 = F( I, J)+S=.M A + S*( B+S'^C) ) 

1 = 1 + 1 

A=(l F( I ,J+1 )-FI I, J ) )+(F( I, J)-F( I ,J-U) )*.5 
0=3.*( F ( I,J+i)-F(I,J))-l2.*A+((F(I,J+2)-F(I,J+i))+(F(I 
1, J + l)-F( I , J) ) )*.5 ) 

C=(A+((F(I,J+2)-F(I,J+l)) + (F(I,J + l)-F(I,J)))=i^.5)~2.=i'(F 
i( I,J + 1)-F( I,J) ) 

F3=F( I , J ) + S-'( A + S*( S+S*C) ) 

1 = 1-3 

A=((F(I,J+1)-F(I,J))+(F(I,J)-F(I,J-1)))*.5 
B = 3.*(Fi I , J+ 1) -F I I , J) ) -(2.* A+( ( F ( I , J+2)-F( I ,J + 1 ) ) + ( F( I 
1, J + l)-F( I, J) ) )-*. 5 ) 

C=(A+((F(I,J+2)-F(I,J+l))+(F(I,J+l)-F(I,J)))«.5)-2.-(F 
1 ( I, J + l )-F( I, J) ) 

F4 = F( I, J ) +S->=( A + S* ( B + S«C ) ) 
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A=((F2-F1)<-(F1-F4) )*.5 

B=3.«{F2-Fl)-(2.'-:'A<-{ (F3-F2) MF2-F1))«.5) 



C=(A+( { F3-F2) + (F2-F1 ) )*. 5 ) -2 .* { F2-F1 ) 
AFF=F1 ( A + R* ( B + R«CJ I 
GO TO 20 

10 AFF=( 1 .-S) =!=( ( 1 .-R )=!=F( I, J) +R*F( H-1 , J ) ) <-S* 
11)4-R*(F( H-1, J<-1) ) ) 

20 RETURN 
END 



( (1 ,-R)=:'F( I, J<- 



SUBROUTINE LLCVT3(GVNI ,GVNJ ,FLAT,FLON) 

IF( (GVNI.GT.30.9.AND.GVM1 ,LT .31 . 1 ) .AND. ( GVN J . GT . 30 . 9 . A 
IND. GVNJ .LT.:>1 . 1 ) ) GO TO 1900 
SETJ=GVNJ-31 . 

SETI=GVNI-31 . 

STEST=0. 

ARTAN=ATAN2( SETJ tSETI ) *57. 295 78 
ATEST=- 10.0 

IF(SETJ .LT.STEST) GO TO 5 

SETK=350.0 

GO TO 10 

5 IF(APTAN.GT.ATEST) GO TO 7 
SETK=ATEST 
GO TO 10 
7 SETK=350.0 
10 FLON=SETK-ARTAN 
SQJ=SET J**2 
SQI=SET 1**2 
R£SQ=973. 752025 

GROUP={ RESQ-SQI-SQJ ) / ( RESO+SQ I + SQ J ) 

FLAT=ARSIN{GRDLIP )*57.2 9578 
GO TO 541 
1000 FLAT=90.0 
FLQN=0. 0 
541 RETURN 
END 



SUBROUTINE H ARM AN ( J 1 , J M, Z , Y , MM ) 

DIME NS I ON FC (36 , 19) , FS( 36, 19) t Y (36 ) , AF( 19) t BF( 19 ) 
IPhASEI 19) ,AMPL( 19) 

IM=36 
II = IM 
IF ( Jl) 10,10,27 
10 J1 = J1 + 1 

Nil = (11/2) +1 

N12 = N1 1 - 1 

El = 6.2831353/FLDAT( ID 

E2 = 1.0/FL0AT(NI12) 

WRITE(6,19) I1,JM,Z,N12 

19 FORMAT! //,' 1 E-W GRID POINTS =',I3,5X, 

I'N-S GRID POINTS =*,I3,5X,' GRID INCREMENT =' 

2, F5. 1 ,5X, 'HIGHEST WAVE NUfMBER POSSIBLE =',I3,//) 



DO 20 I = 1 , IM 
DO 20 J = 1,N11 
FC(I,J) = E2*C3S(FL0AT( ( I 
20 FS(I,J) = E2*SIN( FLOAT! ( I 
DO 25 I = DIM 
FC( I , 1 ) = 0.5*FC ( I , 1 ) 

25 FC(I,N11) = 0.5*FC( I ,N1 1) 
27 N1 1 = (11/2) * 1 
N12 = Nil - 1 
DO 40 J = DNll 
SUM = 0.0 
DO 30 I = DIM 



)*( J-1) )*E1 ) 
)*( J-1 ) )*E1) 



» 
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30 SUM = SUM + Y(I)«FC(I,J) 

40 AF(J) = SUM 

6F(1) = 0.0 
BF(Nll) = 0.0 
DO 60 J = 2,N12 
SUM = 0.0 
DO 50 I = 1, IM 
50 SUM = SUM *■ Y( I )=^FS( I , J) 

60 BF(J) = SUM 

DO 75 J = 1,N11 

PHASE( J ) =0.0 

IF (BF( J) .EO.0.0 ) GO TO 75 

IF ( ( 6FU ) .NE.0.0 ) .AND. ( AF( J) . EQ.0.0) ) GO TO 70 
PHASE(J) = ATAN( 3F(J)/AF( J) ) 

GO TO 75 

70 PHASE(J) = 90.0-(BF( J)/ABS(6F( J) ) ) 

75 CONTINUE 

DO 95 J = IfiNll 

IF (ABS(PHASE( J) ) - 0.8) 80,80,90 
80 AMPL(J) = AF ( J) /COS( PHASEU i ) 

GO TO 95 

90 AMPL(J) = 6F ( J )/SIN( PHASE! J) ) 

95 CONTINUE 

DO 100 J = 1,N11 

100 PHASE(J) = 57.29578*PHASE( J) 

DO 120 J = 2, Nil 
IF (AMPL(J)) 110,120,120 
110 IF (PHASEU) ) 112, lU, 115 

112 PHASE(J) = PHASE(J) + 180 
GO TO 118 

115 PHASE(J) = PHASE(J) - 180 
118 AMPL(J) = -AMPL(J) 

120 CONTINUE 
WRITE(6,121 ) 

121 FORMAT( /,' ',T5,'LAT CIR ‘ , T20, • ARITHMET IC MEAN',T40, 

1 • AMPLITUDE' , T60, • PHASE' ) 

WRITE(6,122) MM,AMPL(1), ((AMPL(I), PH AS E { I ) ) , I =2 , N I 1 ) 

1 22 FORMAT! T6, I2,Tlo,F11.2,5(T36,F8.2,12X,F8.2,/)) 

RETURN 

END 



97 



oooooo ooooooooooooo oonoooo 



ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



THE PURPOSE OF THIS PROGRAM IS TO TAKE A N X N 
MATRIX AND REDUCE IT TO A N X 7 MATRIX TO CONSERVE 
CORE REQUIREMENTS 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 



NPTS= THE GL03AL CORRESPONDENCE TABLE 

NCOR= THE GLOBAL CORRELATION TABLE. THIS IS THE MATRI 
THAT THIS PROGRAM WILL GENERATE. IT CORRELATES THE 
GLOBAL CORRESPCNDENS NUMBER OF THE N X N MATRIX 
WITH IT’S PUSITION IN THE N X 7 MATRIX. 

NTRI5= THE VECTOR THAT CONTAINS THE NUMBER OF THE 
POINT THAT IS SUPPORTED BY ONLY FIVE TRIANGLES VICE 
SIX TRIANGLES. 

NX= A TEMPORY HOLDING VECTOR 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
INTEGER«2 NPTS , NCQ R, NTRI 5 , NX 
COMMON /CMl/ NC3R(642, 7) ,NPTS( 1280,3) 

DIMENSION NTRI 5( 12 ) , I SAVE! 6 ) ,NX ( 7) 

DATA N/ 8/,IK0UNT/l/ 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 



N0TRI= THE NUMBER OF TRIANGLES GIVEN 8 SEGMENTS PER 
MAJOR SPHERICAL TRIANGLE'S SIDE 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
N0TRI=1 280 
DO 1 I=1,NCTRI,5 
L=I+4 

1 READ(5,100)( (NPTS (K,J),J= 1,5) ,K=I,L) 

100 FORMAT! 5( 315) ) 

READ(5,101) NTRI5 

101 FORMAT! 1018) 

N0PTS=642 

DO 20 I=1,N0PTS 
DO 20 J=l,7 
20 NCOR! I , J)=0 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C 

C THE ENTIRE 2 DO LOOP COMPUTES THE GLOBAL CORRELATION 
C TABLE. THE DO 2 LOOP SEARCHES EVERY POINT TO FIND 
C THE FIVE OR SIX TRIANGLE:> THAT SUPPORT THAT POINT. 

C SUBROUTINE CHECK DETERMINES IF THE POINT LIES IN 

C THE TRIANGLE. IF THE POINT LIES IN THE TRIANGLE THEN 

C THE TRIANGLE NUMBER IS STORED IN THE ISAVE VECTOR. 

C THIS IS DONE UNTIL THE PROPER NUMBER OF TRIANGLES 

C HAVE SEEN FOUND. ASSUMING THAT A POINT IS SUPPORTED 

C BY SIX TRIANGLES THE NEXT STEP IS TO SORT OUT THE 

C SEVEN DIFFERENT NUMBERS FROM THE 13 AVAILABLE. 

C SUBROUTINE SORT ACCOMPLISHES THIS AND THEN ARkANGES 

C THE SEVEN NUMBERS IN ASCENDING ORDER. THESE SEVEN 

C NUMBERS ARE THEN PLACED IN THE POINT'S ROW OF THE 

C GLOBAL CORRELATION TABLE. 

C 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DO 2 L=1,N0PTS 
LL = 6 

IF!NTRI5( IKOUNT) .NE.L) GO TO 3 
LL = 5 

IKOU^T=IKOUNT^-l 
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oooooooo 



3 ICCUNT=0 

DO 4 J=1,N0TRI 
IFLAG=0 

CALL CHECK ( NPTS ( J, 1 ) , NPTS( J , 2 ) , N PT S ( J , 3 ) , L , i FL AG ) 

IF( IFLAG.EC.O) GO TO 4 

ICOUNT=lCOUNT<-l 

I SAVE( ICOUNT )=J 

IF ( ICOUNT.EQ.LL ) GO TO 5 

4 CONTINUE 

5 CALL SORT( ISAVE,LL,NX) 

IF(LL.EC.5) LLL=6 
IF(LL.EQ.o) LLL=7 

DO 6 1=1, LLL 

6 NCOR(L, I )=NX( I) 

2 CONTINUE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 



THIS SECTION OUTPUTS THE THE GLOBAL CORRELATION TABLE. 
SINCE THE TABLE NEED BY COMPUTED ONLY ONCE, IT CAN 
BE COMPUTED IN A SEPARATE PROGRAM AND OUTPUT ON 
CARDS AND READ IN AS DATA IN THE MAIN PROGRAM. 



CCCCCCCCCCCCCCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
WPITE(6,102) ( (NCORIK, J) , J=1 ,7) ,K=l,NOPTS) 

102 F0RMAT(/iX,7I10) 

DO 21 I=1,N0PTS 

21 WRITE(7,103)(NC0R( I, J) ,J=1,7) 

103 F0RMATI7I10) 

STOP 

END 



SLBROUT INE CHECK ( OUMl , DUM2 , DUM3 , L , I FLAG ) 

INTEGER*2 OUMl , DUM2 ,DUM3 

IFIDUMl.EQ.U IFLAG=I 

IF(DUM2.EQ.LJ IFLAG=1 

IFIDUM3 .EQ.L ) IFLAG=1 

RETURN 

END 



SUBROUTINE SORT ( I S AVE , LL , NX ) 

INTEGER*2 NPTS , NCO R , NTR 1 5 , NX 

COMMON /CMl/ NCCR(642, 7) ,NPTS( 1280,3) 

DIMENSION NX (7) , I SAVE ( 6) 

NX(1 )=NPTS( ISAVEI 1 ) , 1 ) 

NX(2) =NPTS( I SAVE! 1) ,2) 

NX( 3) =NPTS( I SAVE! 1 ) ,3) 

L = 3 

6 DO 1 1=2, LL 
5 DO 3 K=l,3 
DO 2 J=1,L 

IFINPTSI ISAVEI I) ,K) .EQ.NXIJ ) ) GO TO 3 
IF(J.NE.L) GO TO 2 
NXIL+l ) =NPTS( ISAVEI I ) ,K) 

GO TO 4 

2 CONTINUE 

3 CONTINUE 
GO TO 1 

4 L=L+1 

IFIK.E0.3) GO TO 1 
GO TO 5 
1 CONTINUE 

IFILL.EG.5) LLL=6 
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IF(LL.EQ.6) LLL=7 
DO 8 J=l, 7 
DO 7 I=2 tLLL 

IF(NX( I -1 ) .LT.NXi I ) ) GO TO 7 
ISWAP=NX( I-l ) 

NX( I -1 ) =NX( I J 
NX( I )=ISWAP 

7 CONTINUE 

8 CONTINUE 
RETURN 
END 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



PROGRAM TO GENERATE GEODESIC GRID 
S?RLON= LONGITUDE DF MAJOR SPNtRICAL TRIANGLES 
SPRLAT = LATITUDE OF MAJOR SPHERICAL TRIANGLES 



NPTS(5120,3)= THE GLOBAL CORRESPONDENCE TABLE FOR 
5120 TRIANGLES, N= 16, 

XLAT(6‘t2)= THE LATITUDE VECTOR FOR N= 8 
XL0N1(642) THE LONGITUDE VECTOR FOR N= 8 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
INTEGER*2 NPTS(5120,3) 

DIMENSION SPRLAT ( 12) ,SPRLON( 12) ,XLAT(81 ,49) 

DIMENSION XL0N(31 ,49) 

DIMENSION XLATl ( 642) ,XLQN1( 642 ) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 



THESE FUNCTION STATEMENTS ARE FOR THE LA,^ OF SINES 
(SINLAW), THE LA^ OF COSINES FOR SIDES (COSLAS), AND 
THE LAW OF COSINES FOR ANGLES (COSLAA). ALL THREE 
FUNCTION STATEMENTS ARE FOR SPHERICAL TRIANGLES. 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
S INLAWI A , ANGLA,3 ) = AR S IN ( S I N ( B ) * S I N ( A NG L A ) / S INI A) ) 
COSLASI C, A,ANGLS ) =ARCOS(COS(C)-CCS(A)f 
lSIN(C)=i‘SIN(A)=5‘COS( ANGLB) ) 

COSLAAI ANGLA , ANGLB , ANGLO =ARCOS (ICOS(ANGLB)+ 
1C0S(ANGLC)*C0S( ANGLA) ) 

I/(SIN(ANGLC)=!^SIN( ANGLA) ) ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 



ANGLST= THE INTERIOR ANGLE OF A MAJOR SPHERICAL 
TRIANGLE, 72 DEGREES 
XNOPO= THE LATITUDE OF THE NORTH POLE 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DATA ANGLST/1. 256637062/, N/ 8/ , XNOPO/ 1 . 5 70 79632 7/ 
DATA RAD/57.29577951/ 

NTP1 = 3=!=N+1 

NFF1 = 5*N<-1 

NP1=N+1 

NT = 2=^N 

NF2=N+2 

NM1=N-1 

NM2=N-2 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



NOTRI= NUMBER OF TRIANGLES OVFR THE GLOBE WHERE M= 

THE NUMBER OF SEGMENTS EACH MAJOR SPHERICAL TRIANGLE 
•S SIDES ARE SUBDIVIDED INTO. 

NOPTS= THE NUMBER OF POINTS OVER THE GLOBE 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 





NOTRI=2 


0*>' 


N«=? 


2 




NCPTS=1 


0* 


N«* 


2 + 2 




DO 12 I 


= 1 


,NT 


PI 




DO 12 J 


= 1 


,NF 


PI 




XLONI J , 


I ) 


= 0. 


0 


12 


XLATl J, 


I ) 


= 0. 


0 




READ! 5, 


10 


0) 


SPRLON 




READ! 5, 


10 


0) 


SPRLAT 


100 


FORMAT! 


6F 


12. 


7) 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c 

C SETS UP MODEL PTS. WHERE EACH LONGITUDE LEG IS DIVID- 

C ED INTO N SEGMENTS AND THE PTS. ARE JOINED BY GREAT 

C CIRCES 

C 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

XLON( 1, NPl )=SPRLON(2 ) /RAD 

XLCNI 1,1) =SPKLON ( 1 ) /RAD 

XLCN( NP1,NP1 )=S^RLON(3 )/RAD 

XLAT{ 1, NPl )=SPRLAT(2) /RAD 

XLAT{ 1, 1)=SPRLAT ( 1 )/RAO 

XLAT( NP1,NP1)=S?RI AT(3)/RAD 

ARCLEN=C05LAA( aNGLST, ANGLST , ANGLST) 

SEG=ARCLEN/FLOAT ( N ) 

DO 1 1=2, N 

XLON( 1, I ) =SPRLCN ( 1 ) /RAD 
XLCN{ I , I ) =SPRLON ( 3 ) /RAD 
C0LAT = FL0AT( I-l ) --i'SEG 
XLAT( 1 , I ) = I XNQPO-COLAT) 

1 XLAT< 1 , I ) =XLAT( 1 , I ) 

DC 2 I=2,N 

K=I-1 

ARCLEN=COSLAS(FLOAT( I )«SEG, FLOAT! I )*SEG , ANGLST ) 

ANGLC = S INLAW! ARC LEN, ANGLST, FLOAT! I )-SE3 ) 

ARCL£N = ARCLEN/FLOAT! I ) 

ANGLA=ANGLST/FLOAT( I ) 

DO 2 J=1,K 

XLON! Jf I , H-1 ) =ANGLA'^=FLOAT! J ) 

COLAT = COSLAS!FLOAT! I > =:-'SEG , ARCL E N*FLOAT ! J) , ANGLO 

2 XLAT( Jf 1, I+l )=XNOPO-COLAT 
XLON! l,NTPl)=XLON! 1,1) 

XLAT!1,NTP1)=-XLAT!1,1) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

c 

C THIS DO LOOP MIRRORS THE SOLUTION OF THE FIRST 
C MAJOR SPHERICAL NORTHERN HEMISPHERIC TRIANGLE 
C INTO THE AJOINING FOUR NORTHERN HEMISPHERIC 

C MAJOR SPHERICAL TRIANGLES. 

C 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DO 3 K=l,4 
DO 3 1=2, NPl 
DO 3 J=2,I 

XLAT!K^I+J-K, I )=XLAT( J, I) 

3 XLON! J-K, I ) =XLON( J , I F LOAT ! K ) *ANGL S T 

M=3’:'-N 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

c 

C THIS DO LOOP MIRRORS THE NORTHERN HEMISPHERIC SOLUTION 
C INTO THE SOUTHERN HEMISPERE. 

C 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DO 4 1=2, NPl 
I I=! I-l )«5+l 
DO 14 J = 1 , 1 1 
XLAT! J,M)=-XLAT! J , I ) 

14 XLON! J,M)=XLON! J, I )+- ANGLST/ 2. 

4 M=M-1 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

THIS NEXT SEQUENCE SOLVES THE EQUATORIAL MAJOR 
SPHERICAL TRIANGLES 



ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

A=XN0P0-XL AT ( I,NP1 ) 

C = XNOPO-XLAT( i ,NT+ 1) 

AN&L6=ANGLST/2 . 

ARCLEN=COSLAS(C, A, ANGLB) 

ARCL£N= ARCLEN/ FLOAT! N ) 

M = N 

DO 5 1=1, NMl 

COLAT=COSLAS( A, ARCLcN^FLOAT ( I ) , ANGLST=!'2. ) 

XLAT( 1 , NPI )=XNOPO-COLAT 

XLON( 1,NPUI )=S INLAW (COL AT, ANGLST*2. , ARCLEN=«FLOA T { I ) ) 
XLAT(M,NP1+I)=XLAT(1,NP1+I) 

XLON(M,NPl+I )=ANGLST-XLON( I ,NPH-I) 

5 M=M-1 
M = N 

DC 6 I=1,NM2 
MM1=M-I 

A = XNOPO-XLAT( 1,NP1 + I) 

B=XNGPO-XLAT(M,NPI+I) 

ANGLe=XLON(M,NPl+I )-XLON( 1,NPI+I ) 

ARCLEN=COSLAS ( A , 8 , ANGLB) 

ANGLC=SINLAW( ARC LEN, ANGLB, B) 

IF( I ,LE ,N/2) GO TO 16 
ANGLC= ( XNOPO-ANGLC) +XNOPO 

16 ARCLEN=ARCLEN/FLOAT(MMl ) 

DO 7 J=2,MM1 

COLAT = COSLASi A, ARCLEN^^FLOAT ( J-1 ) , ANGLO 
XLAT( J, MPU-l ) =XNOPO-COLAT 

7 XLOi\'( J,NP1 + I )=SINLAW(COLAT,ANGLC,ARCLEN*FLOAT(J-l) ) + 
1XLCN( l,NPl+I ) 

6 M=N-1 
M = N 

DO 8 1=1, NMl 

XLAT(NPl,NPH-I)=XLAT( 1,NP1 + I) 

8 XLOM( NP1,NP1+I )=XLON{ 1,NP1+I I+ANGLST 
DO 9 I =2, NMl 

ANGLA=XLON(NPl,NPl+I)-XLON( NPl-I ,NP1+I ) 

C = XNOPO-XLAT< NPl - I ,NPH-I ) 

B = XNOPO-XLAT( NPl ,NPH-I ) 

AkCL£N=COSLAS(B, C,ANGLA) 

ANGL6=SINLAW(ARCLEN,ANGLA,C) 

IF( I.LF..N/2) GO TO 17 
ANGLB= ( XNOPO-ANGLB ) +XNOPO 

17 ARCLEN=AkCLEN/FLOAT( I ) 

DO 9 J=2,I 

C0LAT=C0SLAS(B,ARCLEN=:<FL0AT( J-1 ) , ANGLB) 

XLAT(NP2- J,NP1+I )=XNOPO-COLAT 

9 XLON(NP2-J,NPl + I )=XLON(NPl,NPU-I)- 

IS INLAW! COLAT , ANG L 3 , ARCLE N=!=F LOAT ( J-1) ) 

DC 10 K=l,4 
DO 10 1=1, N 
DO 10 J=l,NMl 

XLAT! K=:=N + 1 + I , NPl + J ) =X LAT ! I + 1,NPH-J) 

10 XLCN! K*N+1 + I ,NP1 + J )=XLON! I + 1 ,NPH-J ) ^-FLO A T < K ) ANGL S T 
WRITE!6, 10*^) N,NOPTS 

104 FORMAT! //IX, • THE NUMBER OF PTS IN THE GEODESIC GRID' 
1' DIVIDED INTO' , 2X, 12, 2X, • SEGMENTS IS',2X,I4) 

12X, 12, 2X, ' SEGMENTS IS',2X,I4) 

105 FORMAT! 1X,/24F5. 0) 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



TfiE NEXT SECTION WRITES THE TWO DIMENSIONAL MATRIX 
CONTAINING THE LATITUDES AND THE LONGITUDES INTO A 
ONE DIMENSIONAL VECTOR, ONE FOR THE LATITUDE AND ONE 
FOR THE LONGITUDE. 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
L = 2 

DO 20 J=2,NP1 
MM=( J-1 )*5 
DO 20 1=1, MM 
XLATl ( L )=XLAT( I , J ) 

XLONl ( L)=XLGN( I , J ) 

20 L=L+1 
NTT=NT+1 

DC 41 J=NP2,NTT 
K=5«N 

DO 41 1=1, K 
XLATl ( L) =XLAT( I , J ) 

XLONl ( L)=XLON( I , J) 

L=L+1 
NTT=NTT+1 
KK=NT Pl-1 
DO 21 J=NTT,KK 
K = K-5 

DO 21 I =1 ,K 
XLAT1(L)=XLAT( I, J) 

XLONK L) = XLON( I , J) 

21 L=L+1 

XLATl ( 1 )=XLAT( 1,1) 

XLONK 1 )=XLON( 1,1) 

XLATl (NGPTSJ =XL ATI 1 ,NTP1 ) 

XLONl (NOPTS) =XLON( 1 ,NTP1 ) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 



THIS SECTION DEFINES THE GLOBAL CORRESPONDENCE NUMBER 
FOR EACH POINT IN THE VECTORS CONTAINING THE 
LATITUDES AND THE LONGITUDES. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DO 22 1=1,5 
NPTSI I , 1)=1 
NPTSI I , 2) =1+1 
22 NPTSI I , 3) =1+2 
NPTSI 5,3)=2 
L = 6 
11=2 

ISAVE1=II 
1 1 1=7 
ICCUNT=2 
I SAVE = 7 
DO 23 J=3,NP1 
DO 25 K=l,5 
DO 24 1=1 , ICOUNT 
NPTSI L, 1 )=II 
NPTSIL, 2)=II I 
NFTSIL,3) = II I + l 

I FIK.EQ.5.AND.I.EQ.IC0UNT) GO TO 44 

IFII.EQ. ICOUNT) GO TO 24 

NPTSiL+i , 1 )=I I 

NPTSIL+1,2)=II I+l 

NPTSI L+ 1, 3) = I I + 1 

I FIK.EQ.5 .AND. I . EQ . ICOUNT-1 ) NPTSIL+1,3)=ISAVE1 
I I=I I +1 
I I I = I 1 1 +1 
L = L+2 
GO TO 24 

44 NPTSIL, 1)=ISAVE1 
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■V 



I 



I 



1 



24 

25 

23 



45 

26 

27 



46 

30 

29 



NPTS (L, 3 )=ISAVE 
CCNT INUE 
L = L+1 
I II=I IH-1 
I I = ISAVE 
I SAVE1= I I 
I I I=I SAVE + ( J-1 )«5 
ISAVE = I I I 
ICCUNT= ICOUNT-H 
I II=5*N + I I 
ISAVE=I II 
M=NFP1 -1 
MT=NT+1 

00 27 J=NP2,MT 

00 26 1=1 ,M 
NPTS (L, 1 ) =I I 
NPTS(L,2)=III 
NPTS(L,3) = IH-1 
NPTS(U-1, 1) = I H-1 
NPTSl Lfl , 2 ) = I I I 
NPTS(L<-1 ,3 ) = I I I + l 

I F( I .NS .M) GO TO 45 
NPTS (L, 3)=ISAVE1 
NPTS( U-l, 1 ) = ISAVE1 
NFTS{ U-1 ,3 J= ISAVE 

1 1 = 11 + 1 

I I I = I I I +l 
L = L+2 
CONTINUE 
I I=ISAVe 
I I I = ISAV6 + 5*N 
ISAVE = I I I 
I SAVE1=I I 
MT=MT+1 
NMT=NTP1-1 
ICOUNT= ICOUNT-1 
I SAVE = I II 
KKK=NM1 



00 28 J=MT,MMT 

DO 29 K=l,5 

00 30 1 = 1 , ICOUNT 

NFTSIL, 11=11 

NPTS(L,2)=III 

NPTSa, 3) = I I+l 

I F(K. E0.5 .AND. I . EQ. I COUNT ) 

IFd.EQ. ICOUNT) GO TO 30 

NPTSI L+1 ,11=11+1 

NPTS(L+1,2)=III 

NPTSI L+1, 3 )=I I I+l 

IFIK.EQ.5. AND. I .EQ. ICOUNT- 

I I=I I +1 

I II=I II+l 



L = L+2 



GO TO 50 
NPTS (L, 1 )= ISAVEl 
NPTS(L,2)=ISAVE-1 
CONT INUE 



L = L+1 



11 = 11+1 



I I = ISAVE 
1 I 1 = 1 SAVE + KKK*5 



GO TO 46 



1) NPTS(L+1,3) =ISAVE 



ISAVE1=I I 
KKK=KKK-1 
ISAVE = I I I 

28 IC0UNT=IC0UNT-1 
DC 31 1=1,5 

NPTSINOTR I-5 + 1 , 1 )=N0PTS-6 + I 
NPTS(N0TRI-5 + I , 2 )=NOPTS 
31 NPTS(NuTRI-5+I ,3 )=N0PTS-5+I 
NPTSINOTR I ,3)=NOPTS-5 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c c 
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THE LAST SECTION OUTPUTS THE GLOBAL CORRESPONDENCE 
TABLE, THE LATITUDE AND THE LONGITUDE VECTORS 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
K=N0TRI/6 
DO 32 1=1, K 
I I=K+I 
I I I = 2«K + I 
I I I I = 3’i'Kf I 
I I I II =A«K<-I 
I IIIII = 5=i‘K + I 

32 WRITE(6,110I I , ( N PTS (I , J ) , J = 1 , 3 ) , I I , I NPT S ( I I , J J ) , J J= 1 , 
13),III,(NPTS(III,JJJ),JJJ=1,3),IIII,(NPTS(IIII,JJJJ),J 
2JJJ=1,5),IIIII5(NPTS{IIIII,JJJJJ),JJJJJ=1,3),IIIIII,(N 
3PTS(IIIIII,.JJJJJJ),JJJJJJ = 1,3) 
no FORMAT! /1X,6(^,I5) ) 

WRITE(7,il2) XLATl 
WRITE! 7, I 12) XLONl 
112 FORMAT! 8F10. 6) 

STOP 

END 
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